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Abstract. The aim of the present review is to introduce the reader to some of 
the physical notions and of the mathematical methods that are relevant to the 
study of nonlinear waves in Bose-Einstein Condensates (BECs). Upon introducing 
the general framework, we discuss the prototypical models that are relevant to 
this setting for different dimensions and different potentials confining the atoms. 
We analyze some of the model properties and explore their typical wave solutions 
(plane wave solutions, bright, dark, gap solitons, as well as vortices). We then 
offer a collection of mathematical methods that can be used to understand the 
existence, stability and dynamics of nonlinear waves in such BECs, either directly 
or starting from different types of limits (e.g., the linear or the nonlinear limit, or 
the discrete limit of the corresponding equation). Finally, we consider some special 
topics involving more recent developments, and experimental setups in which there 
is still considerable need for developing mathematical as well as computational 
tools. 
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Abbreviations 

• AC: Anti-Continuum (Limit) 

• BEC: Bosc-Einstcin condensate 

• BdG: Bogoliubov-de Gennes (Equations) 
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• cqNLS: cubic-quintic NLS 

• DNLS: Discrete Nonlinear Schrodinger (Equation) 

• EP: Ermakov-Pinncy (Equation) 

• GP: Gross-Pitaevskii (Equation) 

• KdV: Korteweg-de Vries (Equation) 

• LS: Lyapunov-Schmidt (Technique) 

• MT: Magnetic Trap 

• NLS: Nonlinear Schrodinger (Equation) 

• NPSE: Non-polynomial Schrodinger Equation 

• ODE: Ordinary Differential Equation 

• OL: Optical Lattice 

• PDE: Partial Differential Equation 

• RPM: Reductive Perturbation Method 

• TF: Thomas-Fermi 

1. Introduction 

The phenomenon of Bose-Einstein condensation is a quantum phase transition 
originally predicted by Bose [1] and Einstein [2,3] in 1924. In particular, it was 
shown that below a critical transition temperature T c , a finite fraction of particles 
of a boson gas (i.e., whose particles obey the Bose statistics) condenses into the 
same quantum state, known as the Bose-Einstein condensate (BEC). Although Bose- 
Einstein condensation is known to be a fundamental phenomenon, connected, e.g., 
to superfluidity in liquid helium and superconductivity in metals (see, e.g., Ref. [4]), 
BECs were experimentally realized 70 years after their theoretical prediction: this 
major achievement took place in 1995, when different species of dilute alkali vapors 
confined in a magnetic trap (MT) were cooled down to extremely low temperatures 
[5-7], and has already been recognized through the 2001 Nobel prize in Physics [8,9]. 
This first unambiguous manifestation of a macroscopic quantum state in a many- 
body system sparked an explosion of activity, as reflected by the publication of several 
thousand papers related to BECs since then. Nowadays there exist more than fifty 
experimental BEC groups around the world, while an enormous amount of theoretical 
work has followed and driven the experimental efforts, with an impressive impact on 
many branches of Physics. 

From a theoretical standpoint, and for experimentally relevant conditions, the 
static and dynamical properties of a BEC can be described by means of an effective 
mean-field equation known as the Gross-Pitaevskii (GP) equation [10,11]. This is a 
variant of the famous nonlinear Schrodinger (NLS) equation [12] (incorporating an 
external potential used to confine the condensate), which is known to be a universal 
model describing the evolution of complex field envelopes in nonlinear dispersive media 
[13]. As such, the NLS equation is a key model appearing in a variety of physical 
contexts, ranging from optics [14-17], to fluid dynamics and plasma physics [18], 
while it has also attracted much interest from a mathematical viewpoint [12,19,20]. 
The relevance and importance of the NLS model is not limited to the case of 
conservative systems and the theory of solitons [13,18,21-23]; in fact, the NLS equation 
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is directly connected to dissipative universal models, such as the complex Ginzburg- 
Landau equation [24,25], which have been studied extensively in the context of pattern 
formation [26] (see also Ref. [27] for further discussion and applications). 

In the case of BECs, the nonlinearity in the GP (NLS) model is introduced by the 
interatomic interactions, accounted for through an effective mean-field. Importantly, 
the mean-field approach, and the study of the GP equation, allows the prediction and 
description of important, and experimentally relevant, nonlinear effects and nonlinear 
waves, such as solitons and vortices. These, so-called, matter-wave solitons and 
vortices can be viewed as fundamental nonlinear excitations of BECs, and as such have 
attracted considerable attention. Importantly, they have also been observed in many 
elegant experiments using various relevant techniques. These include, among others, 
phase engineering of the condensates in order to create vortices [28,29] or dark matter- 
wave solitons in them [30-34], the stirring (or rotation) of the condensates providing 
angular momentum creating vortices [35,36] and vortex-lattices [37-39], the change of 
scattering length (from repulsive to attractive via Feshbach resonances) to produce 
bright matter-wave solitons and soliton trains [40-43] in attractive condensates, or 
set into motion a repulsive BEC trapped in a periodic optical potential referred to as 
optical lattice to create gap matter- wave solitons [44]. As far as vortices and vortex 
lattices are concerned, it should be noted that their description and connection to 
phenomena as rich and profound as superconductivity and superfluidity, were one of 
the themes of the Nobel prize in Physics in 2003. 

The aim of this paper is to give an overview of some physical and mathematical 
aspects of the theory of BECs. The fact that there exist already a relatively large 
number of reviews [45-51] and textbooks [4,52-55] devoted in the Physics of BECs, 
and given the space limitations of this article, will not allow us to be all-inclusive. 
Thus, this review naturally entails a personalized perspective on BECs, with a special 
emphasis on the nonlinear waves that arise in them. In particular, our aim here is to 
present an overview of both the physical setting and, perhaps more importantly, of the 
mathematical techniques from dynamical systems and nonlinear dynamics that can be 
used to address the dynamics of nonlinear waves in such a setting. This manuscript 
is organized as follows. 

Section 2 is devoted to the mean-field description of BECs, the GP model and its 
properties. In particular, we present the GP equation and discuss its variants in the 
cases of repulsive and attractive interatomic interactions and how to control them via 
Feshbach resonances. We also describe the ground state properties of BECs and their 
small-amplitude excitations via the Bogoliubov-de Gennes equations. Additionally, 
we present the types of the external confining potential and how their form leads to 
specific types of simplified mean-field descriptions. 

Section 3 describes the reduction of the spatial dimensionality of the BEC by 
means of effectively suppressing one or two transverse directions. This can be achieved 
by "tightening" the external confining potential (usually a harmonic magnetic trap) 
along these directions. We introduce the basic nonlinear structures (dark and bright 
solitons) that are ubiquitous to one-dimensional settings. The different types of 
nonlinearities that arise from different approximations due to the dimensionality 
reduction are discussed. We also present the dimensionality reduction in the presence 
of external periodic potentials generated by the optical lattices (which are created 
as interference patterns of multiple laser beams) and the discrete limit, the discrete 
nonlinear Schrodinger equation, that they entail for strong potentials. 

Section 4 deals with the mathematical methods used to describe nonlinear waves 
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in BECs. The presentation concerns four categories of methods, depending on the 
particular features of the model at hand. The first one concern "direct" methods, 
which analyze the nonlinear mean-field models directly, without employing techniques 
based on some appropriate, physically relevant and mathematically tractable limit. 
Such approaches include, for example the method of moments, self-similarity and 
rescaling methods, or the variational techniques among others. The second one will 
concern methods that make detailed use of the understanding of the linear limit of 
the problem (e.g., the linear Schrodinger equation in the presence of a parabolic, 
periodic, or a double- well potential). The third category of the mathematical methods 
entails perturbation techniques from the fully nonlinear limit of the system (e.g., the 
intcgrable NLS equation), while the fourth one concerns discrete systems (relevant 
to BECs trapped in strong optical lattices), where perturbation methods from the 
so-called anti-continuum limit are extremely helpful. 

Finally, in Sec. 5 we present some special topics that have recently attracted much 
physical interest, both theoretical and experimental. These include multicomponent 
and spinor condensates described by systems of coupled GP equations, shock waves, as 
well as nonlinear structures arising in higher-dimensions, such as vortices and vortex 
lattices in BECs, and multidimensional solitons (including dark and bright ones). We 
also briefly discuss the manipulation of matter-waves by means of various techniques 
based on the appropriate control of the external potentials. In that same context, 
the effect of disorder on the matter-waves is studied. Finally, we touch upon the 
description of BECs beyond mean-field theory, presenting relevant theoretical models 
that have recently attracted attention. 



2. The Gross-Pitaevskii (GP) mean-field model 

2.1. Origin and fundamental properties of the GP equation 

We consider a sufficiently dilute ultracold atomic gas, composed by N interacting 
bosons of mass m confined by an external potential V cy ±{r). Then, the many-body 
Hamiltonian of the system is expressed, in second quantization form, through the 
boson annihilation and creation field operators, ^(r,*) and W{v,t), as [46,53], 

H = y"dr* t (r,t)H *(r,t) + ^ J drdr'tftfr, t)&{r', f)V(r - r')*(r', t)tf(r, t), (1) 

where Ho = —(h 2 /2m)V 2 + T4 xt (r) is the "single-particle" operator and V(r — r') 
is the two-body interatomic potential. The mean-field approach is based on the 
so-called Bogoliubov approximation, first formulated by Bogoliubov in 1947 [56], 
according to which the condensate contribution is separated from the boson field 
operator as $>(r,t) = ^(r, t) + w(r,t). In this expression, the complex function 
$?(r,t) = (^(r, i)) (i.e., the expectation value of the field operator), is commonly 
known as the macroscopic wavefunction of the condensate, while w(r',t) describes 
the non-condensate part, which, at temperatures well below T c , is actually negligible 
(for generalizations accounting for finite temperature effects see Sec. 5.7). Then, the 
above prescription leads to a nontrivial zeroth-order theory for the BEC wavefunction 
as follows: First, from the Heisenberg evolution equation ih(d^> /dt) = [*&,-/?] for the 
field operator (r, t) , the following equation is obtained: 



i^*(r,f) 



H Q + J rfr'* t (r',t)F(r'-r)*(r',t) 



*(r,t). (2) 



CONTENTS 



6 



Next, considering the case of a dilute ultracold gas with binary collisions at low energy, 
characterized by the s-wave scattering length a, the interatomic potential can be 
replaced by an effective delta-function interaction potential, V(r' — r) = g5(r' — r) 
[46,47,52,53], with the coupling constant (i.e., the nonlinear coefficient) g given by 
g = 4irh 2 a/m. Finally, employing this effective interaction potential, and replacing 
the field operator ^ with the classical field 'J, Eq. (2) yields the GP equation, 



- — V 2 + K xt (r)+. 9 |v]/(r,i)| : 



tt(r,t). 



(3) 



The complex function * in the GP Eq. (3) can be expressed in terms of the 
density p(r,t) = |<J/(r,£)| 2 , and phase S(r,t) of the condensate as \I>(r, t) = 
y/p(r,t) exp[iS(r, t)}. Note that the current density j = 2^j(**V*-*V**) (asterisk 
denotes complex conjugate), assumes a hydrodynamic form j = pv, with an atomic 
velocity v(r, t) = —VS(r,t). The latter is irrotational (i.e., Vxv = 0), which is a 
typical feature of superfluids, and satisfies the famous Onsager-Feynman quantization 
condition § c dl ■ v = (h/m)Af, where Af is the number of vortices enclosed by the 
contour C (the circulation is obviously zero for a simply connected geometry). 

For time-independent external potentials, the GP model possesses two integrals 
of motion, namely, the total number of atoms, 



N = j |*(r,t)| 2 dr, 



and the energy of the system, 



E 



-I 



dv 



I 



2m 



W\ z + V c ^\<S>\ z + -g\y 



(4) 



(5) 



with the three terms in the right-hand side representing, respectively, the kinetic 
energy, the potential energy and the interaction energy. 

A time-independent version of the GP Eq. (3) can readily be obtained upon 
expressing the condensate wave function as ^>(r,t) = ^o(r) exp(— ipst/H), where 
is a function normalized to the number of atoms (N — J dr l^ol 2 ) and fi = dE/dN 
is the chemical potential. Substitution of the above expression into the GP Eq. (3) 
yields the following steady state equation for ^ ■ 



V 2 

2m 



+ K x t(r)+ 5 |*o| 2 (r) 



*o(r) = /i*o(r). 



(6) 



Equation (6) is useful for the derivation of stationary solutions of the GP equation, 
including the ground state of the system (see Sec. 2.5). 



2.2. The GP equation vs. the full many-body quantum mechanical problem 

It is clear that the above mean-field approach and the analysis of the pertinent GP 
Eq. (3) is much simpler than a treatment of the full many-body Schrodinger equation. 
However, a quite important question is if the GP equation can be derived rigorously 
from a self-consistent treatment of the respective many-body quantum mechanical 
problem. Although the GP equation is known from the early 1960s, this problem 
was successfully addressed only recently for the stationary GP Eq. (6) in Ref. [57] . In 
particular, in that work it was proved that the GP energy functional describes correctly 
the energy and the particle density of a trapped Bose gas to the leading-order in the 



CONTENTS 



7 



small parameter p|a| 3 , where p is the average density of the gas.^f The above results 
were proved in the limit where the number of particles N — > oo and the scattering 
length a — ► 0, such that Na is fixed. Importantly, although Ref. [57] referred to the 
full three-dimensional (3D) Bose gas, extensions of this work for lower-dimensional 
settings were also reported (see the review [51] and references therein). 

The starting point of the analysis of Ref. [57] is the effective Hamiltonian of TV 
identical bosons. Choosing the units so that ft = 2m = 1, this Hamiltonian is expressed 
as (see also Ref. [52]), 

N 

i=i i<3 

where u(|r|) is a general interaction potential assumed to be spherically symmetric and 
decaying faster than |r|~ 3 at infinity. Then, denoting the quantum-mechanical ground- 
state energy of the Hamiltonian (7) (which depends on the number of particles N and 
the dimensionless + scattering length a) by Eqm(N,o), the main theorem proved in 
Ref. [57] is as follows: 

• The GP energy is the dilute limit of the quantum-mechanical energy: 

V5i > : lim ±-E QM (n, —) = Egp(1,oi), (8) 

where £"gp {N, a) is the energy of a solution of the dimcnsionlcss stationary GP 
Eq. (6) (in units such that ft = 2m = 1), and the convergence is uniform on 
bounded intervals of a\. 

The above results (as well as the ones in Ref. [51]) were proved for stationary 
solutions of the GP equation, and, in particular, for the ground state solution. More 
recently, the time-dependent GP Eq. (3) was also analyzed within a similar asymptotic 
limit (N — > oo) in Ref. [58]. In this work, it was proved that the limit points of the k- 
particlc density matrices of *jv,t (which is the solution of the TV-particle Schrodingcr 
equation) satisfy asymptotically the GP equation (and the associated hierarchy of 
equations) with a coupling constant given by J v(x)dx, where v(x) describes the 
interaction potential. 

Thus, these recent rigorous results justify (under certain conditions) the use of the 
mean-field approach and the GP equation as a quite relevant model for the description 
of the static and dynamic properties of BECs. 



2.3. Repulsive and attractive interatomic interactions. Feshbach resonance 

Depending on the BEC species, the scattering length a [and, thus, the nonlincarity 
coefficient g in the GP Eq. (3)] may take either positive or negative values, accounting 
for repulsive or attractive interactions between the atoms, respectively. Examples of 
repulsive (attractive) BECs are formed by atomic vapors of 87 Rb and 23 Na ( 85 Rb and 
7 Li) ones, which are therefore described by a GP model with a defocusing (focusing) 
nonlincarity in the language of nonlinear optics [12,17]. 

% When JV|a| 3 <C 1, the Bose-gas is called "dilute" or "weakly-interacting". In fact, the smallncss 
of this dimcnsionlcss parameter is required for the derivation of the GP Eq. (3); in typical BEC 
experiments this parameter takes values iV|a| 3 < 10 — 3 [46]. 

+ The dimensionless two-body scattering length is obtained from the solution u(r) of the zero- 
energy scattering equation —u"(r) + ^v(r)u(r) = with u(0) = and is given, by definition, as 
a = lim (r — u(r)/u'(r)) (see also Refs. [52,53]). 
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On the other hand, it is important to note that during atomic collisions, the 
atoms can stick together and form bound states in the form of molecules. If the 
magnetic moment of the molecular state is different from the one of atoms, one may 
use an external (magnetic, optical or dc-electric) field to controllably vary the energy 
difference between the atomic and molecular states. Then, at a so-called Feshbach 
resonance (see, e.g., Ref. [59] for a review), the energy of the molecular state becomes 
equal to the one of the colliding atoms and, as a result, long-lived molecular states are 
formed. This way, as the aforementioned external field is varied through the Feshbach 
resonance, the scattering length is significantly increased, changes sign, and finally 
is decreased. Thus, Feshbach resonance is a quite effective mechanism that can be 
used to manipulate the interatomic interaction (i.e., the magnitude and sign of the 
scattering length). 

Specifically, the behavior of the scattering length near a Feshbach resonant 
magnetic field B is typically of the form [60,61], 



where a is the value of the scattering length far from resonance and A represents 
the width of the resonance. Feshbach resonances were studied in a series of 
elegant experiments performed with sodium [62,63] and rubidium [64,65] condensates. 
Additionally, they have been used in many important experimental investigations, 
including, among others, the formation of bright matter- wave solitons [40-43]. 

2.4- The external potential in the GP model 

The external potential V eKt (r) in the GP Eq. (3) is used to trap and/or manipulate the 
condensate. In the first experiments, the BECs were confined by means of magnetic 
fields [8,9], while later experiments demonstrated that an optical confinement of BECs 
is also possible [66,67], utilizing the so-called optical dipole traps [68,69]. While 
magnetic traps are typically harmonic (see below), the shape of optical dipole traps 
is extremely flexible and controllable, as the dipole potential is directly proportional 
to the light intensity field [69]. An important example is the special case of periodic 
optical potentials called optical lattices (OLs), which have been used to reveal novel 
physical phenomena in BECs [70-72]. 

In the case of the "traditional" magnetic trap, the external potential has the 
harmonic form: 



where, in general, the trap frequencies u x , oj y , oj z along the three directions are 
different. On the other hand, the optical lattice is generated by a pair of laser beams 
forming a standing wave which induces a periodic potential. For example, a single 
periodic ID standing wave of the form E(z, t) — 2E cos(fcz) exp(— iujt) can be created 
by the superposition of the two identical beams, E± (z,t) = E n exp[i(±kz— cot)] , having 
the same polarization, amplitude E 0} wavelength A = 2n/k 7 and frequency ui. Since 
the dipole potential Vdi P is proportional to the intensity I ~ \E(z, t)\ 2 of the light field 
[69], this leads to an optical lattice of the form Vdi P = Vol = Vb cos 2 (fcz). In such a 
case, the lattice periodicity is A/2 and the lattice height is given by Vb <~ / max /Aa;, 
where I max is the maximum intensity of the light field and Aui = uj — uj is the detuning 
of the lasers from the atomic transition frequency ui . Note that atoms are trapped 




(9) 




m(uj x x 2 + uj y y 2 + uj z z 2 ), 



(10) 
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at the nodes (anti-nodes) of the optical lattice for blue- (red-) detuned laser beams, 
or Auj > (Aw < 0). In a more general 3D setting, the optical lattice potential can 
take the following form: 

VolM = V [cos 2 (k x x + <p x ) + cos 2 (k y y + <f> y ) + cos 2 (k z z + <p z )] , (11) 

where h = 27r/Ai (i e {x,y, z}), Aj = A/[2sin(0j/2)], 9i are the (potentially variable) 
angles between the laser beams [70,72] and are arbitrary phases. 

It is also possible to realize experimentally an "optical superlattice" , characterized 
by two different periods. In particular, as demonstrated in Ref. [73], such a superlattice 
can be formed by the sequential creation of two lattice structures using four laser 
beams. A stationary ID superlattice can be described as V(z) — V\ cos(k\z) + 
Vi cos(fc2z), where fej and Vi denote, respectively, the wavenumbers and amplitudes of 
the sublattices. The experimental tunability of these parameters provides precise and 
flexible control over the shape and time-variation of the external potential. 

The magnetic or/and the optical dipole traps can be experimentally combined 
cither together, or with other potentials; an example concerns far off-resonant laser 
beams, that can create effective repulsive or attractive localized potentials, for blue- 
detuned or red-detuned lasers, respectively. Such a combination of a harmonic trap 
with a repulsive localized potential located at the center of the harmonic trap is the 
double-well potential, as, e.g., the one used in the seminal interference experiment 
of Ref. [74]. Double- well potentials have also been created by a combination of 
a harmonic and a periodic optical potential [75]. Finally, other combinations, 
including, e.g., linear ramps of (gravitational) potential V cx t — mgz have also been 
experimentally applied (see, e.g., Refs. [75,76]). Additional recent possibilities include 
the design and implementation of external potentials, offered, e.g., by the so-called 
atom chips [77-79] (see also the review [80]). Importantly, the major flexibility for the 
creation of a wide variety of shapes and types of external potentials (e.g., stationary, 
time-dependent, etc), has inspired many interesting applications (see, for example, 
Sec. 5.4). 

2.5. Ground state 

The ground state of the GP model of Eq. (3) can readily be found upon expressing the 
condensate wave function as \I/(r,i) = ^>o(r) exp(—i/j,t/h). If g = 0, Eq. (6) reduces 
to the usual Schrodinger equation with potential V cx t- Then, for a harmonic external 
trapping potential [see Eq. (10)], the ground state of the system is obtained when 
letting all non-interacting bosons occupy the lowest single-particle state; there, <I>o 
has the Gaussian profile 



where Who = (w x w v uj z ) 1 /' i is the geometric mean of the confining frequencies. 

For repulsive interatomic forces (g > 0, or scattering length a > 0), if the number 
of atoms of the condensate is sufficiently large so that Na/aho ^> 1, the atoms are 
pushed towards the rims of the condensate, resulting in slow spatial variations of the 
density. Then the kinetic energy (gradient) term is small compared to the interaction 
and potential energies and becomes significant only close to the boundaries. Thus, 
the Laplacian kinetic energy term in Eq. (6) can safely be neglected. This results in 



3/4 



/ 2 2 i 2 2 , 2 2\ 

cxp --Ki +uj y y +u z z ) 



(12) 
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the, so-called, Thomas-Fermi (TF) approximation [46,53,52] for the system's ground 
state density profile: 

p(r) = |* (r)| 2 =g- 1 [»-V ext (r)}, (13) 

in the region where \x > V ext (r), and p = outside. For a spherically symmetric 
harmonic magnetic trap (V cxt = Vmt with ui^o = u) x = w y = w z ), the radius 
-Rtf = (2/i/?n) 1/ ' 2 /a;ho for which p(Rtf) = 0, is the so-called Thomas-Fermi 
radius determining the size of the condensed cloud. Furthermore, the normalization 
condition connects fj, and N through the equation /j, = (hw ho /2)(15Na/a^ ) 2 ^ 5 , where 
iho = (fr/m-^ho) 1 ^ 2 is the harmonic oscillator length. 

For attractive interatomic forces (g < 0, or a < 0), the density tends to increase 
at the trap center, while the kinetic energy tends to balance this increase. However, if 
the number of atoms N in the condensate exceeds a critical value, i.e., N > N cr , the 
system is subject to collapse in 2D or 3D settings [12,52,53]. Collapse was observed 
experimentally in both cases of the attractive 7 Li [81] and 85 Rb condensate [82]. In 
these experiments, it was demonstrated that during collapse the density grows and, as 
a result, the rate of collisions (both elastic and inelastic) is increased; these collisions 
cause atoms to be ejected from the condensate in an energetic explosion, which leads to 
a loss of mass that results in a smaller condensate. It should be noted that the behavior 
of BECs close to collapse can be significantly affected by effects such as inelastic two- 
and three-body collisions that are not included in the original GP equation; such 
effects are briefly discussed below (see Sec. 3.4). 

The critical number of atoms necessary for collapse in a spherical BEC is 
determined by the equation iV cr |a|/aho = 0.575 [83], where \a\ is the absolute value of 
the scattering length. Importantly, collapse may not occur in a quasi-lD setting (see 
Sec. 3.1), provided that the number of atoms does not exceed the critical value given 
by the equation N cr \a\/a r = 0.676, with a r being the transverse harmonic oscillator 
length [84-86]. 



2.6. Small- amplitude linear excitations 



We now consider small-amplitude excitations of the condensate, which can be 
found upon linearizing the time-dependent GP equation around the ground state. 
Specifically, solutions of Eq. (3) can be sought in the form 



*(r,t) = e 



*o(r) + J2 (%(r)e" <W3 '* + v*{v)e^ 1 ) 



(14) 



where Uj, Vj are small (generally complex) perturbations, describing the components 
of the condensate's (linear) response to the external perturbations that oscillate at 
frequencies ±u)j [the latter are (generally complex) eigenfrequencies] . Substituting 
Eq. (14) into Eq. (3), and keeping only the linear terms in Uj and vj, the following 
set of equations is derived 

H Q - i u + 2 3 |*o| 2 (r) Uj {r) + g * 2 (r)^-(r) = hu jUj (r), 

Ho - n + 2g |* | 2 (r)j w,-(r) + g% 2 (r) Uj (r) - -huj Vj (r), (15) 

where Hq = —(h 2 /2m)V 2 + V ex t(r). Equations (15) are known as the Bogoliubov-de 
Gennes (BdG) equations. These equations can also be derived using a purely quantum- 
mechanical approach [46,47,52,53] and can be used, apart from the ground state, 
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for any state (including solitons and vortices) with the function being modified 
accordingly. 

The BdG equations are intimately connected to the stability of the state \I>o- 
Specifically, suitable combinations of Eqs. (15) yield 

(ivj-cj*) j(\u ] \ 2 -\v, J \ 2 )dv = Q. (16) 

This equation can be satisfied in two different ways: First, if uj — u* — 0, i.e., if the 
eigenfrequencies ujj are real; if this is true for all j, the fact that Im{cjj} = shows 
that the state is stable. Note that, in this case, one can use the normalization 
condition for the cigcnmodes Uj, Vj of the form J{\uj\ 2 — \vj\ 2 )dr = 1. On the other 
hand, occurrence of imaginary or complex eigenfrequencies uj (i.e., if uj — uo* ^ or 
Imjcijj} 7^ 0), indicates dynamical instability of the state <J/o> m such a case, Eq. (16) 
is satisfied only if J \uj\ 2 dv = J \vj\ 2 dr. 

In the case of a uniform gas (i.e, for V cxt (r) = and = P =const.), the 
amplitudes u and v are plane waves <~ e jk r (of wavevector k) and the BdG Eqs. (15) 
lead to a dispersion relation, known as the Bogoliubov spectrum: 

w =(^j(^r + H- (17) 

For small momenta hk, Eq. (17) yields the phonon dispersion relation w — cq, where 
c=\fgpjm (18) 

is the speed of sound, while, for large momenta, the spectrum provides the free particle 
energy h 2 k 2 / (2m); the "crossover" between the two regimes occurs when the excitation 
wavelength is of the order of the healing length [see Eq. (19)]. 

In the case of attractive interatomic interactions (g < 0), the speed of sound 
becomes imaginary, which indicates that long wavelength perturbations grow or decay 
exponentially in time. This effect is directly connected to the modulational instability, 
which leads to delocalization in momentum space and, in turn, to localization in 
position space and the formation of solitary-wave structures. Modulational instability 
is responsible for the formation of bright matter- wave solitons [40-42] , as was analyzed 
by various theoretical works (see, e.g., Refs. [87-89] and the reviews [90,91]). 



3. Lower-dimensional BECs, solitons, and reduced mean-field models 

3.1. The shape of the condensate and length scales 

In the case of the harmonic trapping potential (10), the flexibility over the choice of 
the three confining frequencies uij (j € {x, y, z}) may be used to control the shape of 
the condensate: if ui x = uj v = ui r w uj z the trap is isotropic and the BEC is almost 
spherical, while the cases lo z < uj r or uj r < uj z describe anisotropic traps in which the 
BEC is, respectively, "cigar shaped" , or "disk-shaped" . The strongly anisotropic cases 
with lo z <§C uj r or uj r <§C uj z are particularly interesting as they are related to effectively 
quasi one-dimensional (ID) and quasi two-dimensional (2D) BECs, respectively. Such 
lower dimensional BECs have been studied theoretically [92-98] (see also Ref. [51] for 
a rigorous mathematical analysis) and have been realized experimentally in optical 
and magnetic traps [99], in optical lattice potentials [100-103] and surface microtraps 
[78,79]. 
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The confining frequencies of the harmonic trapping potential set characteristic 
length scales for the spatial size of the condensate through the characteristic harmonic 
oscillator lengths aj = {h/mojj) 1 ^ 2 . Another important length scale, introduced by the 
effective mean-field nonlinearity, is the healing length, which is the distance over which 
the kinetic energy and the interaction energy balance: if the BEC density grows from 
to p over the distance £, the kinetic energy, <~ h 2 /(2m^ 2 ), and interaction energy, 
~ 4:Trh 2 ap/m, become equal at the value of £ given by [46,52,53] 

f= (8Trpa)- 1/2 . (19) 

Note that the name of £ is coined by the fact that it is actually the distance over which 
the BEC wavefunction "heals" over defects. Thus, the spatial widths of nonlinear 
excitations, such as dark solitons and vortices in BECs, arc of 0(£). 



3.2. Lower-dimensional GP equations 

Let us assume that io z <C co x = u> y = u> r . Then, if the transverse harmonic oscillator 
length a r = ^h/muu r < £, the transverse confinement of the condensate is so tight 
that the dynamics of such a cigar-shaped BEC can be considered to be effectively ID. 
This allows for a reduction of the fully 3D GP equation to an effectively ID GP model, 
which can be done for sufficiently small trapping frequency ratios uo z /uo r - It should be 
stressed, however, that such a reduction should be only considered as the ID limit of a 
3D mean-field theory and not as a genuine ID theory (see, e.g., Ref. [51] for a rigorous 
mathematical discussion). Similarly, a disk-shaped BEC with a z < £ and sufficiently 
small frequency ratios u r /ui z , can be described by an effective 2D GP model. Below, 
we will focus on cigar-shaped BECs and briefly discuss the case of disk-shaped ones. 

Following Refs. [84,104,105] (see also Ref. [91]), we assume a quasi-lD setting 
with lo z <C uv and decompose the wavefunction ^ in a longitudinal (along z) and a 
transverse [on the (x,y) plane] component; then, we seek for solutions of Eq. (3) in 
the form 

tf(r,i)=#M)$(r;t), (20) 

where $(r;i) = $o(r) ex P(~ *7^)> r2 = % 2 +V 2 i while the chemical potential 7 and the 
transverse wavefunction $(r) are involved in the auxiliary problem for the transverse 
quantum harmonic oscillator, 

h 2 ~ 1 

— V*$ - -mu 2 r rH + 7*o - 0, (21) 

where V 2 = d 2 /dx 2 + d 2 /dy 2 . Since the considered system is effectively ID, it is 
natural to assume that the transverse condensate wavefunction $(r) remains in the 
ground state; in such a case <&o( r ) takes the form $>o(r) = n^^-^a^ 1 cxp(— r 2 j2a 2 r ) 
[note that when considering the reduction from 3D to 2D the transverse wave function 

takes the form $ (r) = 7r" 1 / 4 a7 1/2 exp(-r 2 /2a^)]- Then, substituting Eq. (20) into 
Eq. (3) and averaging the resulting equation in the r-direction (i.e., multiplying by 
and integrating with respect to r), we finally obtain the following ID GP equation, 



ifr^ip{z,t) 



■~+V{z) +giD Mz,t)\> 



1>{z,t), (22) 



where the effective ID coupling constant is given by gm = g/2na^. 
V(z) — (l/2)mco 2 z 2 . On the other hand, in the 2D case of disk-shaped BECs, 
the respective (2 + l)-dimcnsional NLS equation has the form of Eq. (22), with d 2 
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being replaced by the Laplacian V 2 ., the effectively 2D coupling constant is g 2 D = 
g/V2na z — 2\/27raa z huj z , while the potential is V(x,y) = (l/2)muj 2 (x 2 + y 2 ). Note 
that such dimensionality reductions based on the averaging method are commonly 
used in other disciplines, as, e.g., in nonlinear fiber optics [17]. 

A similar reduction can be performed if, additionally, an optical lattice potential is 
present. In this case, it is possible (as, e.g., in the experiment of Ref. [44]) to tune lo z so 
that it provides only a very weak trapping along the z-direction; this way, the shift in 
the potential trapping energies over the wells where the BEC is confined can be made 
practically negligible. In such a case, the potential in Eq. (22) is simply the ID optical 
lattice V(z) = Vo cos 2 (kz). Similarly, in the quasi-2D case, an "egg-carton potential" 
V(x,y) = Vq [cos 2 (k x x) + cos 2 (k y y)] is relevant for disk-shaped condensates. 

We note in passing that the dimensionality reduction of the GP equation can also 
be done self-consistently, using multiscale expansion techniques [106,107]. It is also 
worth mentioning that more recently a rigorous derivation of the ID GP equation 
was presented in Ref. [108], using energy and Strichartz estimates, as well as two 
anisotropic Sobolev inequalities. 



3.3. Bright and dark matter-wave solitons 

The ID GP Eq. (22) can be reduced to the following dimensionless form, 



2 



-^ + V(z)+g\iP(z,t)\ 2 



i/>(z,t), (23) 



where the density |f/'| 2 , length, time and energy are respectively measured in units 
of 47r|a|a 2 ., a r , uj^ 1 and fou) r , while the coupling constant g is rescaled to unity (i.e., 
g = ±1 for repulsive and attractive interatomic interactions respectively). In the case 
of a homogeneous BEC (V(z) = 0), Eq. (23) becomes the "traditional" completely 
intcgrable NLS equation. The latter, is well-known (see, e.g., Ref. [21]) to possess an 
infinite number of conserved quantities (integrals of motion), with the lowest-order 
ones being the number of particles: 



the momentum: 



and the energy: 



N= / 



p = (i/2) (tpr z -r^),dz 



/ — oo 
(hM 2 + sM 4 ) dz, 
-oo 



where the subscripts denote partial derivatives. 

The type of soliton solutions of the NLS equation depends on the parameter g. In 
particular, for attractive BECs (g = —1), the NLS equation possesses a bright soliton 
solution of the following form [109], 

ipbs(z,t) = t]sech[r](z — vt)]exp[i(kz — ut)], (24) 

where 77 is the amplitude and inverse spatial width of the soliton, while k, to and 
v = duj/dk — k are the soliton wavenumber, frequency, and velocity, respectively. 
The frequency and wavenumber of the soliton are connected through the "soliton 
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dispersion relation" uj — i(fc 2 — rj 2 ), which implies that the allowable region in the 
(k,ui) plane for bright solitons is located below the parabola u> = \k 2 1 corresponding 
to the "elementary excitations" (i.e., the linear wave solutions) of the NLS equation. 

Introducing the solution (24) into the integrals of motion N, P and E it is readily 
found that 

N = 2r,, P = 2r,k, E = r,k 2 - \r, 3 . (25) 

These equations imply that the bright soliton behaves as a classical particle with 
effective mass M\, s , momentum Pbs and energy E^s, respectively given by Mb s = 2r], 
-Pbs = MbsV, and Eb s — ■^M\ >s v 2 — -^M^ s , where it is reminded that v — k. Notice 
that in the equation for the energy, the first and second terms in the right hand 
side are, respectively, the kinetic energy and the binding energy of the quasi-particles 
associated with the soliton [110]. Differentiating the soliton energy and momentum 
over the soliton velocity, the following relation is found, 

m: = v > (26) 

which underscores the particle-like nature of the bright soliton. 

On the other hand, for repulsive BECs (g = +1), the NLS equation admits 
a dark soliton solution, which in this case lives on the nonzero background ip = 
ipoexp[i(kz — ujt)]. The dark soliton may be expressed as [111], 

ip(z,t) — ipo (cos ip tanh ( + i sin ip) exp[i(kz — u)t)] , (27) 

where ( = ip cosip i z — vt), v = (l/2)^ 2 + ipo, while the remaining parameters v, ip 
and k, are connected through the relation v = tpo sm( P + k. Here, tp is the so-called 
"soliton phase angle", or, simply, the phase shift of the dark soliton (\ip\ < tt/2), which 
describes the darkness of the soliton through the relation, \^\ 2 = 1 — cos 2 <^sech 2 £; this 
way, the limiting cases <p = and cosy <C 1 correspond to the so-called black and 
gray solitons, respectively. The amplitude and velocity of the dark soliton are given 
by cosy and siny respectively; thus, the black soliton, tp = ipo tanh(^oa;) exp(— i/it), 
is a stationary dark soliton (v = 0), while the gray soliton moves with a velocity close 
to the speed of sound (v <~ c = ipo in our units). The dark soliton solution (27) has 
two independent parameters, for the background (ipo and k) and one for the soliton 
(ip). In fact, it should be mentioned that in both the bright and the dark soliton, 
there is also a freedom in selecting the initial location of the solitary wave z (in the 
above formulas, zq has been set equal to zero) *. Also, it should be noted that as in 
this case the dispersion relation implies that lo > k , the allowable region in the (k, u>) 
plane for dark solitons is located above the parabola w = \k 2 . 

As the integrals of motion of the NLS equation refer to both the background and 
the dark soliton, the integrals of motion for the dark soliton are renormalized so as to 
extract the contribution of the background [112-114]. In particular, the renormalized 
momentum and energy of the dark soliton (27) read (for k = 0): 

{c 2_ v 2y/2- 



P ds = - 2v(c 2 - v 2 ) 1/2 + 2c 2 tan 



(28) 



E ds = A -{c 2 -v 2 f/ 2 . (29) 

* Recall that the underlying model, namely the completely integrable NLS equation, has infinitely 
many symmetries, including translational and Galilean invariances. 
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Upon differentiating the above expressions over the soliton velocity v, it can readily 
be found that 

dE ds 

mr v > m 

which shows that, similarly to the bright soliton, the dark soliton effectively behaves 
like a classical particle. Note that, usually, dark matter-wave solitons are considered 
in the simpler case where the background is at rest, i.e., k — 0; then, the frequency 
u actually plays the role of a normalized one-dimensional chemical potential, namely 
/i = tpQ, which is determined by the number of atoms of the condensate. Moreover, 
it should be mentioned that in the case of a harmonically confined condensate, 
i.e., for V(z) = \^z 2 (with ft = uo z /uo r being the normalized trap strength) in 
Eq. (23), the background of the dark soliton is actually the ground state of the BEC 
which can be approximated by the Thomas-Fermi cloud [see Eq. (13)]; thus, the 
"composite" wavefunction (containing both the background and the dark soliton) can 
be approximated e.g. by the form ip = iPtf(z) exp(— iip 2 t)ipds(z , t), where Vds(z,i) is 
the dark soliton of Eq. (27). 

Both types of matter-wave solitons, namely the bright and the dark ones, have 
been observed in a series of experiments. In particular, the formation of quasi-lD 
bright solitons and bright soliton trains has been observed in 7 Li [40,41] and 85 Rb [42] 
atoms upon tuning the interatomic interaction within the stable BEC from repulsive 
to attractive via the Feshbach resonance mechanism (discussed in Sec. 3.3). On the 
other hand, quasi-lD dark solitons were observed in 23 Na [30,31] and 87 Rb [32-34] 
atoms upon employing quantum-phase engineering techniques or by dragging a moving 
impurity (namely a laser beam) through the condensate. Note that multidimensional 
solitons (and vortices), as well as many interesting applications based on the particle- 
like nature of matter-wave solitons highlighted here will be discussed in Sec. 5. 



3.4- Mean-field models with non-cubic nonlinearities 

Mean-field models with non-cubic nonlinearities have also been derived and used in 
various studies, concerning either the effect of dimensionality on the dynamics of 
cigar-shaped BECs, or the effect of the three-body collisions irrespectively of the 
dimensionality of the system. 

Let us first discuss the former case, i.e., consider a condensate confined in a highly 
anisotropic trap with, e.g., lo z -C u> r and examine the effect of the deviation from 
one-dimensionality on the longitudinal condensate dynamics. Following Refs. [115— 
119], one may factorize the wavefunction as per Eq. (20), but with the transverse 
wavefunction $ depending also on the longitudinal variable z. Then, one may employ 
an adiabatic approximation to separate the fast transverse and slow longitudinal 
dynamics (i.e., neglecting derivatives of $ with respect to the slow variables z and t). 
This way, assuming that the external potential is separable, Kxt( r ) = U(r) + V(z), 
the following system of equations is obtained from the 3D GP Eq. (3), 

ih tt=-^U + V ^ + ^ ^ 

/v(p)$ = -|^V2$ + tf( r )$ + .gp|$| 2 $, (32) 
2m 

where the transverse local chemical potential fJ, r {p) (which depends on the longitudinal 
density p(z, t) — \i()(z, t)\ 2 ) is determined by the ground state solution of Eq. (32). An 
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approximate solution of the above system of Eqs. (31)-(32) was found in Ref. [115] 
(see also Refs. [116,117]) as follows. As the system is close to ID, it is natural to 
assume that the transverse wave function is close to the ground state of the transverse 
harmonic oscillator, and can be expanded in terms of the radial eigenmodes $j, i.e., 
<E>(r; z) — &o(r)+J2j Cj(z)$j(r). Accordingly, expanding the chemical potential /v(p) 
in terms of the density as /x r (p) = hw r + gip — g 2 p 2 + • • •, the following NLS equation 
is obtained: 



dtp 



1>, (33) 



with the nonlinearity function given by 

f(p)=9ip-92P 2 - (34) 

It is clear that Eq. (33) is a cubic-quintic NLS (cqNLS) equation, with the coefficient 
of the linear and quadratic term being given by Ref. [115] g\ = gm = 2ahuj r and 
#2 = 24 ln(4/3)a 2 ftuv, respectively. In the effectively ID case discussed in Sec. 3.2, 
this cqNLS equation is reduced to the ID GP Eq. (22). Note that the cqNLS model has 
been used in studies of the dynamics of dark [115] and bright [116,117] matter-wave 
solitons in elongated BECs. 

The transverse chemical potential p r of an elongated condensate was also derived 
recently by Muhoz Mateo and Delgado [118] as a function of the longitudinal density 
p. This way, the same authors presented in the recent work of Ref. [119] the effective 
ID NLS Eq. (33), but with the nonlinearity function given by 



f(p) = y/i + AaNp, (35) 

where a is the scattering length and N is the number of atoms. jj Note that in the 
weakly- interacting limit, 4aNp <C 1, the resulting NLS equation has the form of 
Eq. (22), with the same coupling constant gm- This model, which was originally 
suggested in Ref. [120], predicts accurately ground state properties of the condensate, 
such as the chemical potential, the axial density profile and the speed of sound. 

Other approaches to the derivation of effective lower-dimensional mean-field 
models have also been proposed in earlier works, leading (as in the case of 
Refs. [118-120]) to NLS-type equations with generalized nonlinearities [121-125]. 
Among these models, the so-called non-polynomial Schrodinger equation (NPSE) 
has attracted considerable attention. The latter was obtained by Salasnich et al. 
[121] by employing the following ansatz for the transverse wavefunction, <f>(r; t) = 
[exp (— r 2 /2a 2 {z, £))] /[TT 1 ^ 2 a(z, t)]; then, the variational equations related to the 
minimization of the action functional (from which the 3D GP equation can be derived 
as the associated Euler-Lagrange equation) led to Eq. (33) with a nonlinearity function 

and to the following equation for the transverse width a: a 2 = a 2 ^/l + 2aNp. Note 
that in the weakly interacting limit of aNp <C 1, Eq. (33) becomes again equivalent 
to the ID GP Eq. (22), while the width a becomes equal to the transverse harmonic 
oscillator length a r . The NPSE (33) has been found to predict accurately static and 
dynamic properties of cigar-shaped BECs (such as the density profiles, the speed 



U Note that the number of atoms N appears in the nonlinearity function f(p) due to the fact that 
the wavefunction is now normalized to 1 rather than to N, as in the GP Eq. (22). 



CONTENTS 



17 



of sound, and the collapse threshold of attractive BECs) [126], while its solitonic 
solutions have been derived in Ref. [121]. Generalizations of the NPSE model in 
applications involving time-dependent potentials [124,127], or the description of spin-1 
atomic condensates [125], have also been presented. Moreover, it is worth mentioning 
that the NPSE has been found to predict accurately the BEC dynamics in recent 
experiments [75]. 

On the other hand, as mentioned above, mean-field models with non-cubic 
nonlinearities, and particularly the cqNLS equation in a ID, 2D or a 3D setting, 
may have a different physical interpretation, namely to take into account three-body 
interactions. In this context, and in the most general case, the coefficients g\ and g 2 
in Eq. (34) are complex, with the imaginary parts describing inelastic two- and three- 
body collisions, respectively [128]. As concerns the three-body collision process, it 
occurs at interparticle distances of order of the characteristic radius of interaction 
between atoms and, generally, results in the decrease of the density that can be 
achieved in traps. Particularly, the rate of this process is given by (dp/dt) = —Lp 3 
[52], where p is the density and L is the loss rate, which is of order of 10 -27 -10~ 30 
cm 6 s _1 for various species of alkali atoms [129]. Accordingly, the decrease of the 
density is equivalent to the term — (L/2)|\f f | 4 \t r in the time dependent GP equation, 
i.e., to the above mentioned quintic term. 

The cqNLS model has been studied in various works, mainly in the context 
of attractive BECs (scattering length a < 0, or g\ < 0). Particularly, in studies 
concerning collapse, both cases with real g 2 [130] , and with imaginary g 2 [131,132] were 
considered. Additionally, relevant lower-dimensional (and in particular ID) models 
were also analyzed in Refs. [133,134]. The latter works present also realistic values 
for these ID cubic-quintic NLS models, including also estimations for the three-body 
collision parameter g 2 (see also Rcfs. [135-137] in which the coefficient g 2 was assumed 
to be real). Additionally, periodic potentials were also considered in such cubic- 
quintic models and various properties and excitations of the BECs, such as ground 
state and localized excitations [135], or band-gap structure and stability [136], were 
studied. Moreover, studies of modulational instability in the continuous model with 
the dissipative quintic term [138], or the respective discrete model with a conservative 
quintic term [137] were also reported. 

3.5. Weakly- and strongly-interacting ID Bose gases. The Tonks- Girardeau gas 

In the previous subsections we discussed the case of ultracold weakly-interacting quasi- 
1D BECs, which, in the absence of thermal or quantum fluctuations, are described 
by an effectively ID GP equation [cf. Eq. (22)]. On the other hand, and in the same 
context of ID Bose gases, there exists the opposite limit of strong interatomic coupling 
[92,93,139]. In this case, the collisional properties of the bosonic atoms are significantly 
modified, with the interacting bosonic gas behaving like a system of free fermions; such, 
so-called, Tonks- Girardeau gases of impenetrable bosons [140,141] have recently been 
observed experimentally [142,143]. The transition between the weakly and strongly 
interacting regimes is usually characterized by a single parameter 7 — 2/ (pain) [92,93] , 
with aiD = Or/ a 3-D an d 03.0 being the effective ID and the usual 3D scattering lengths, 
respectively (a r is the transverse harmonic oscillator length) [53]. This parameter 
quantifies the ratio of the average interaction energy to the kinetic energy calculated 
with mean-field theory. Notice that 7 varies smoothly as the interatomic coupling is 
increased from values 7 <C 1 for a weakly interacting ID Bose gas, to 7 » 1 for the 
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strongly-interacting Tonks-Girardeau gas, with an approximate "crossover regime" 
being around 7 ~ O(l), attained experimentally as well [144,145]. 

The above mentioned weakly- and strongly-interacting ID Bose gases can 
effectively be described by a generalized ID NLS of the form of Eq. (33). In 
such a case, while the functional dependence of f{p) on 7 (and its analytical 
asymptotics) are known [146], its precise values in the crossover regime can only be 
evaluated numerically. Such intermediate values have been tabulated in Ref. [139], 
and subsequently discussed by various authors in the framework of the local density 
approximation [147,148]. Note that following the methodology of Ref. [121], Salasnich 
et al. have proposed a model different from the NPSE, but still of the form of Eq. (33), 
with the nonlinearity function depending on the density |^| 2 and the transverse width 
a of the gas, to describe the weakly- and the strongly-interacting regimes, as well as 
the crossover regime [149]. Finally, it is noted that a much simpler approximate model 
(with f(p) being an explicit function of the density), which also refers to these three 
regimes, was recently proposed in Ref. [150]. 

Coming back to the case of the Tonks-Girardeau gas, it has been suggested that 
an effective mean-field description of this limiting case can be based on a ID quintic 
NLS equation, i.e., an equation of the form (33) with a nonlinearity function [151]: 

/w -^-.. . . (37) 

The quintic NLS equation was originally derived by Kolomcisky et al. [152] from a 
renormalization group approach, and then by other groups, using different techniques 
[153-155]; it is also worth noticing that its time-independent version has been 
rigorously derived from the many-body Schrodinger equation [156]. Although the 
applicability of the quintic NLS equation has been criticized (as in certain regimes 
it fails to predict correctly the coherence properties of the strongly-interacting ID 
Bose gases [157]), the corresponding hydrodynamic equations for the density p and 
the phase S arising from this equation are well-documented in the context of the local 
density approximation [139,147]. In fact, this equation is expected to be valid as 
long as the number of atoms exceeds a certain minimum value (typically much larger 
than 10), for which oscillations in the density profiles become essentially suppressed 
[151,154,158]; in other words, the density variations should occur on a length scale 
which is larger than the Fermi healing length £p = l/(np p ) (where p p is the peak 
density of the trapped gas). 

The quintic NLS model has been used in various studies [151,159-161], basically 
connected to the dynamics of dark solitons in the Tonks-Girardeau gas, and in the 
aforementioned crossover regime of 7 ~ 0(1) [150]. In this connection, it is relevant 
to note that in the above works it was found that, towards the strongly- interacting 
regime, the dark soliton oscillation frequency is up-shifted, which may be used as a 
possible diagnostic tool of the system being in a particular interaction regime. 



3. 6. Reduced mean-field models for BECs in optical lattices 

Useful reduced mean-field models can also be derived in the case where the BECs are 
confined in periodic (optical lattice) potentials. Here, we will discuss both continuous 
and discrete variants of such models, focusing — as in the previous subsections — on 
the ID case (generalizations to higher-dimensional settings will be discussed as well). 
We start our exposition upon considering that the external potential in Eq. (22) is 
of the form V(z) — Vosin 2 (/cz), i.e., an optical lattice of periodicity L = ir/k. Then, 
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measuring length, energy and time in units of a L = L/ir, E L = 2E rcc = ti 2 /ma 2 L 
(where the recoil energy E rcc is the kinetic energy gained by an atom when it absorbs 
a photon from the optical lattice), and w^ 1 = fr/E L , respectively, we express Eq. (22) 
in the dimensionless form of Eq. (23) with V(z) = sin 2 (z). 

Generally, the stationary states of Eq. (23) can be found upon employing the usual 
ansatz, ^(z,t) — F(z) exp(— i^t), where \x is the dimensionless chemical potential. 
In the limiting case g — > (i.e., for a noninteracting condensate) where the Bloch- 
Floquet theory is relevant, the function F(x) can be expressed as [162] F(z) = 
Uk,a(z) exp(ikz), where the functions Uk, a {z) share the periodicity of the optical 
lattice, i.e., Uk, a {z) — Uk^ a {z + nL) where n is an integer. If the Floquet exponent 
(also called "quasimomcntum" in the physics context) k is real, the wavefunction ip 
has the form of an infinitely extended wave, known as a Block wave. Such waves exist 
in bands (which are labeled by the index a introduced above), while they do not exist 
in gaps, which are spectral regions characterized by Im(/c) ^ 0. 

The concept of Bloch waves can also be extended in the nonlinear case (g ^ 0) 
[163-170]. In particular, when the coupling constant is small, the nonlinear band- 
gap spectrum and the nonlinear Bloch waves are similar to the ones in the linear 
case [166]. However, when the coupling constant is increased (or, physically speaking, 
the local BEC density grows), the chemical potential of the nonlinear Bloch wave is 
increased (decreased) for g > (g < 0) and, thus, the nonlinearity effectively "shifts" 
the edges of the linear band. In this respect, it is relevant to note that for a sufficiently 
strong nonlinearity, "swallowtails" (or loops) appear in the band-gap structure, both 
at the boundary of the Brillouin zone and at the zone center. This was effectively 
explained in Ref. [170], where an adiabatic tuning of a second lattice with half period 
was considered. Swallowtails in the spectrum are related to several interesting effects, 
such as a non-zero Landau-Zener tunneling probability [163], the existence of two 
nonlinear complex Bloch waves (which are complex conjugate of each other) at the 
edge of the Brillouin zone [165,166] , as well as the existence of period-doubled states (in 
the case of sufficiently strong optical lattices - see Sec. 3.6.2), also related to periodic 
trains of solitons [169] in this setting. We finally note that experimentally it is possible 
to load a BEC into the ground or excited Bloch state with an unprecedented control 
over both the lattice and the atoms [103]. 



3.6.1. Weak optical lattices. Let us first consider the case of weak optical lattices, 
with Vq <C yU- I n this case, and in connection to the above discussion, a quite relevant 
issue is the possibility of nonlinear localization of matter- waves in the gaps of the linear 
spectrum, i.e., the formation of fundamental nonlinear structures in the form of gap 
solitons, as observed in the experiment [44] . The underlying mechanism can effectively 
be described by means of the so-called Bloch-wave envelope approximation near the 
band edge, first formulated in the context of optics [171], and then used in the BEC 
context as well [172,106,173] (see also Ref. [174] and the review [175]). In particular, 
a multiscale asymptotic method was used to show that the BEC wave function can 
effectively be described as ip(z,t) = U{z,t)ukfi{z)cyi^>{ikz), where Ukfi{z) exp(ikz) 
represents the Bloch state in the lowest band a = (at the corresponding central 
quasimomentum fc), while the envelope function U{z,t) is governed by the following 
dimensionless NLS equation: 



i- 



.dU _ 1 d 2 U 
dt 2m e g dz 2 



9'id\U\ 2 U, (38) 
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where g' 1D is a renormalized (due to the presence of the optical lattice) effectively 
fD coupling constant, and m c g is the effective mass. Importantly, the latter is 
proportional to the inverse effective diffraction coefficient d 2 /j,/dk 2 whose sign may 
change, as it is actually determined by the curvature of the band structure of the linear 
Bloch waves. Obviously, the NLS Eq. (38) directly highlights the abovementioned 
nonlinear localization and soliton formation, which occurs for m C Hg' 1D < 0. Note 
that the above envelope can be extended in higher-dimensional (2D and 3D) settings 
[176], in which the effective mass becomes a tensor. In such a case, and for repulsive 
BECs (g < 0), all components of the tensor have to be negative for the formation of 
multi-dimensional gap solitons [177,178]. 

Coupled-mode theory, originally used in the optics context [171,179], has also been 
used to describe BECs in optical lattices [180-184]. According to this approach, and 
in the same case of weak optical lattice strengths, the wavefunction is decomposed into 
forward and backward propagating waves, A(z, t) and B(z, t), with momenta k = +1 
and k = — 1, respectively, namely 

ip(z, t) — [A(z, t) exp(ix) + B(z, t) exp(— ix)] exp(— i/it). 

This way, Eq. (23) can be reduced to the following system of coupled-mode equations 
(see also the recent work [185] for a rigorous derivation), 

l (^ + 2 ^) =VoB + g(\Af + 2\Bf)A, (39) 

1 (lF ~ 2 S) = VoA + 9(21 A? + m 

Notice that Eqs. (39)-(40) are valid at the edge of the first Brillouin zone, i.e., in the 
first spectral gap of the underlying linear problem (with g — 0). There, the assumption 
of weak localization of the wave function (which is written as a superposition of just 
two momentum components) is quite relevant: for example, a gap soliton near the 
edge of a gap can indeed be approximated by a modulated Bloch wave, which itself is 
a superposition of a forward and backward propagating waves [171]. Thus, coupled- 
mode theory was successfully used to describe matter-wave gap solitons in Refs. [180- 
184]. Note that the coupled-mode Eqs. (39)-(40) can directly be connected by a NLS 
equation of the form (38) by means of a formal asymptotic method [183] that uses 
the distance from the band edges as a small parameter. We finally mention that the 
coupled-mode equations can formally be extended in higher dimensions [185]. 

3.6.2. Strong optical lattices and the discrete nonlinear Schrodinger equation. 
Another useful reduction, which is relevant to deep optical lattice potentials with 
Vo 3> ^, is the one of the GP equation to a genuinely discrete model, the so- 
called discrete NLS (DNLS) equation [186]. Such a reduction has been introduced 
in the context of arrays of BEC droplets confined in the wells of an optical lattice in 
Refs. [187,188] and further elaborated in Ref. [189]; we will follow the latter below. 

When the optical lattice is very deep, the strongly spatially localized wave 
functions at the wells of the optical lattice can be approximated by Wannicr 
functions, i.e., the Fourier transform of the Bloch functions. Due to the completeness 
of the Wannier basis, any solution of Eq. (23) can be expressed as ip(z,t) = 
J2 n a c n.a{t)w nta (z), where n and a label wells and bands, respectively. Substituting 
the above expression into Eq. (23), and using the orthonormality of the Wannier basis, 
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we obtain a set of differential equations for the coefficients. Upon suitable decay of the 
Fourier coefficients and the Wannier functions' prefactors (which can be systematically 
checked for given potential parameters), the model can be reduced to 
dc 

— w 0,« c «,« + (C n -l,a + C n +l,a) 

+ 9 W / Q 0l a 2 03 C n,a 1 c ",t>2 C n,t>3! (41) 

ai ,a2 ,ft3 

where W^l^lal = J- 00 Wn,a'U)n 1 ,a 1 'U) n2>a2 w n3 , a3 dx. The latter equation degenerates 
into the so-called tight-binding model [187,188], 

dc n a 

i—lf- = <^o, a c n ,a + ^i,a {c n -i, a + c n+lia ) + gW^l n \c nia \ 2 c nia , (42) 

if one restricts consideration only to the first band. Equation (42) is precisely the 
reduction of the GP equation to its discrete counterpart, namely, the DNLS equation. 
Higher-dimensional versions of the latter are of course physically relevant models 
and have, therefore, been used in various studies concerning quasi-2D and 3D BECs 
confined in strong optical lattices (see subsequent sections). 



4. Some mathematical tools for the analysis of BECs 

Our aim in this section will be to present an overview of the wide array of mathematical 
techniques that have emerged in the study of BECs. Rather naturally, one can envision 
multiple possible partitions of the relevant methods, e.g., based on the type of the 
nonlinearity (repulsive or attractive, depending on the sign of the scattering length), or 
based on the type of the external potential (periodic, decaying or confining) involved 
in the problem. However, in the present review, we will classify the mathematical 
methods based on the mathematical nature of the underlying considerations. We 
will focus, in particular, on four categories of methods. The first one will concern 
"direct" methods which analyze the mean-field model directly, without initiating 
the analysis at some appropriate, mathematically tractable limit. Such approaches 
include, e.g., the method of moments, self-similarity and rescaling methods, or the 
variational techniques among others. The second one will concern methods that make 
detailed use of the understanding of the linear limit of the problem (in a parabolic, or a 
periodic potential or in combinations thereof including, e.g., a double- well potential). 
The third will entail perturbations from the nonlinear limit of the system (such as, 
e.g., the integrable NLS equation), while the fourth one will concern discrete systems 
where perturbation methods from the so-called anti-continuum limit of uncoupled sites 
are extremely helpful. 



4-1. Direct methods 

Perhaps one of the most commonly used direct methods in BEC is the so-called 
variational approximation (see Ref. [190] for a detailed review). It consists of using an 
appropriate ansatz, often a solitonic one or a Gaussian one (for reasons of tractability 
of the ensuing integrations) in the Lagrangian or the Hamiltonian of the model at hand, 
with some temporally dependent variational parameters. Often these parameters 
are the amplitude and/or the width of the BEC wavefunction. Then, subsequent 
derivation of the Euler-Lagrange equations leads to ordinary differential equations 
(ODEs) for such quantities which can be studied either analytically or numerically 
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shedding light on the detailed dynamics of the BEC system. Such methods have 
been extensively used in examining very diverse features of BECs including collective 
excitations [191], studying the dynamics of BECs in ID optical lattices [187], offering 
insights on the collapse or absence thereof in higher dimensional settings and potentials 
— see, e.g., Ref. [192] and references therein — , or on the behavior of BECs on space- or 
time-dependent nonlinearity settings — see e.g., Ref. [193] as an example — . However, 
both due to the very widespread use of the method and due to the fact that detailed 
reviews of it already exist [190,194], we will not focus on reviewing it here. Instead, 
we direct the interested reader to the above works and references therein. We also 
note in passing one of the concerns about the validity of the variational method, which 
consists of its strong restriction of the infinite-dimensional GP dynamics to a small 
finite dimensional subspace (freezing the remaining directions by virtue of the selected 
ansatz). This restriction is well-known to potentially lead to invalid results [195]; it is 
worthwhile to note, however, that there are efforts underway to systematically compute 
corrections to the variational approximation [196], thereby increasing the accuracy of 
the method. 

Another very useful tool for analyzing BEC dynamics is the so-called moment 
method [197], whereby appropriate moments of the wavefunction ijj — ^/p exp(i<fi) 
(where p = |i/>| 2 and <p are the BEC density and phase, respectively) are defined such 
as N — J pdr (the number of atoms), Xi = J Xipdr (the center of mass location), 
Vi = J pd(p/dxidr (the center speed), Wi = J x 2 pdr (the width of the wavefunction), 
Bi = 2 j Xipd(f>/dxidr (the growth speed), K t = —(1/2) J ip*d 2 4>/dx 2 dr (the kinetic 
energy) and J = J G(p) dr (the interaction energy). Notice that in the above 
expressions the subscript i denotes the i-th direction. Then for the rather general 
GP-type mean-field model of the form: 



with a generalized nonlinearity g{p) — dG/dp, the generalized dissipation a and, 
say, the typical parabolic potential of the form V(r) = J2k^/^) U} t x k^ tnc mom ent 
equations read [197]: 




^AV + V(r)V + g(\ip\ 2 ,t)iP - ta(\TP\ 2 ,t)iP, 



(43) 




(44) 




(45) 




(46) 




(47) 




(48) 




(49) 
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dJ ^ f ^a 2 (t> , n f . f dG 



(50) 



where SG = G(p)—pg(p). One can then extract, for parabolic potentials, a closed-form 
exact ODE describing the motion of the center of mass (assuming that the dissipation 
does not dependent on p) of the form: 

±^L = - L o 2 (t)X l -2a(t)^-2a(t)X l . (51) 
dt z dt 

In the absence of dissipation (and for constant in time magnetic trap frequencies), 

this yields a simple harmonic oscillator for the center of mass of the condensate. This 

is the so-called Kohn mode [198], which has been observed experimentally (see, e.g., 

Refs. [52,53]). In fact, more generally, for conservative potentials one obtains a general 

Newtonian equation of the form [199]: 

= - / i^rP dr ( 52 ) 



dt 2 J dxi 

which is the analog of the linear quantum-mechanical Ehrcnfest theorem. 

There are some simple dissipationless (i.e., with a = 0) cases for which the 
Eqs. (44)-(50) close. For example, if the potential is spherically symmetric (u>i(t) = 
Lo(t)), one can close the equations for R = VW, together with the equation for K into 
the so-called Ermakov-Pinney (EP) equation [200] of the form: 

d 2 R M 

= -«W* +13. «*) 

where M is a constant depending only on the initial data and the interaction strength 
U (the equations close only for the two-dimensional case and for G = Up 2 ). One of the 
remarkable features of such EP equations [200] is that they can be solved analytically 
provided that the underlying linear Schrodinger problem d 2 R/dt 2 = —u(t)R is 
explicitly solvable with linearly independent solutions Ri(t) and R2(t). In that case, 
the EP equation has a general solution of the form (AR 2 + BR?, + 2CR1R2) 1 / 2 , where 
the constants satisfy AB — C 2 — Mj w 2 and w is the Wronskian of the solutions Ri 
and i?2- Such EP equations have also been used to examine the presence of parametric 
resonances for time-dependent frequencies (such that the linear Schrodinger problem 
has parametric resonances) in Refs. [201,202]. Another place where such EP approach 
has been quite relevant is in the examination of BECs with temporal variation of the 
scattering length; the role of the latter in preventing collapse has been studied in the 
EP framework in Ref. [203] . It has also been studied in the context of producing exact 
solutions of the second moment of the wavefunction, which is associated with the 
width of the BEC, which are either oscillatory (breathing condensates) or decreasing 
in time (collapsing condensates) or increasing in time (dispersing BECs) in Ref. [204] . 
It should be mentioned that when the scattering length is time-dependent the EP 
equation (53) is no longer exact, but rather an approximate equation, relying on the 
assumption of a quadratic spatial dependence of the condensate phase. Another case 
where exact results can be obtained for the moment equations is when the nonlincarity 
g is time-independent and the phase satisfies Laplace's equation Acf> = 0, in which case 
vortex-line solutions can be found for the wavefunction ip [197]. 

It should also be mentioned in passing that such moment methods are also used 
in deriving rigorous conditions for avoiding collapse in NLS equations more generally, 
where this class of methods is known under the general frame of variance identities 
(see, e.g., the detailed discussion of Sec. 5.1 in Ref. [12]). 
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Another comment regarding the above discussion is that Newtonian dynamical 
equations of the form of Eq. (52) are more generally desirable in describing the 
dynamics not only of the full wavefunction but also of localized modes (nonlinear 
waves), such as bright solitons. This approach can be rigorously developed for small 
potentials V(x) — eW(x) or wide potentials V(x) — W(ex) in comparison to the 
length scale of the soliton. In such settings, it can be proved [205] that the motion of 
the soliton is governed by an equation of the form 

= - VC/ 00> (54) 

where m e s is an effective mass (found to be 1/2 independently of dimension in 
Ref. [205]) and the effective potential is given by 

In ID, Eq. (55) can be used to characterize not only the dynamics of the solitary 
wave, but also its stability around stationary points such that U'(so) = 0. It is 
natural to expect that its motion will comprise of stable oscillations if U"(s n ) > 0, 
while it will be unstable for U"(sq) < 0. In the case of multiple such fixed points, the 
equation provides global information on the stability of each equilibrium configuration 
and local dynamics in the neighborhood of all equilibria. In higher dimensions, an 
approach such as the one yielding Eq. (54) is not applicable due to the instability of 
the corresponding multi-dimensional (bright) solitary waves to collapse [12]. This type 
of dynamical equations was originally developed formally using asymptotic multi-scale 
expansions as, e.g., in Refs. [206,207]. We will return to a more detailed discussion of 
such techniques characterizing the dynamics of the nonlinear wave in Sees. 4.3.2 and 
4.3.3. 

Another class of methods that can be used to obtain reduced ODE 
information from the original GP partial differential equation (PDE) concerns scaling 
transformations, such as the so-called lens transformation [12]. An example of this 
sort with 

4>(x, t) = J^ cxp(z/(t)r 2 ) u r(t)j , (56) 

has been used in Ref. [208] to convert the more general (with time dependent 
coefficients) form of the GP equation 

^ = -^T A ^ + 5 fi (*) r V + 9W ~ (57) 

into the simpler form with time independent coefficients: 

.du 1 1|2 < . 

i— — = — -Ajjit + s\u\ u, (58) 

where n = x/L{t) and s = ±1. This happens if the temporally dependent functions 
l(t), f(t), L(t) and r(t) satisfy the similarity conditions: 

^ = a(t)dl + a(t)l, (59) 
| = -2a{t)f-\m, (60) 
§ =Mt).fL, (61) 
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dr a(t) 



dt L 2 



(62) 



= ■ (63) 

Some of these ODEs can be immediately solved, e.g., 

L(t)= exp^*a(0/w) , (64) 

i(t) = r(i)exp^d jf a(t')f(t')dt'^ =T{t)L(t)i, (65) 

ff (t) = S a(t)r(t)L(t) d - 2 , (66) 

where T(t) — exp(/ * a(t')dt'). Notice that this indicates that a(t),g(t) and a(t) are 
inter-dependent through Eq. (66). While, unfortunately, Eq. (60) cannot be solved in 
general, it can be integrated in special cases, such as, e.g., Q(i) = 0. Notice that in 
that case, periodic a(t) with zero average, i.e., a(t')dt' — implies that L(t),l(t) 
and f(t) will also be periodic. In other settings the above equations can be used to 
construct collapsing solutions or to produce pulsating two-dimensional profiles, as is 
the case, e.g., for Q(t) = m(l — 2sn 2 (t, m)) in the form: 

1 ( x \ 2 mcn(t, m)sn(i, m) \ , . 

* = M^f { dn(t,m),r(t) ) ^ " 2dnft, to) J ' (6?) 

where U(rf) is the well-known 2D Townes soliton (see also Sec. 5.4.2) profile and 
r(t) — f Q dn(t', m)~ 2 dt'. Similar types of lens transformations were used to study 
collapse type phenomena [209,210] and to examine modulational instabilities in the 
presence of parabolic potentials [87]. 

In the same as the above class of transformation methods one can classify also the 
scaling methods that arise from the consideration of Lie group theory and canonical 
transformations [211] of nonlinear Schrodingcr equations with spatially inhomogencous 
nonlinearities of the form (for the stationary problem) 

-ip xx + V(x)^+g(x)^ 3 = M V- (68) 

Considering the generator of translational invariance motivates the scaling of the form 
U(x) = bix)- 1 ' 2 ^ and X = j*[l/b(s)}ds with g(x) = g /b(x) 3 . Then, U satisfies 
the regular ID NLS equation (whose solutions are known from the inverse scattering 
method [21]) and b satisfies: 

b"'(x) - 2b(x)V'(x) + 4b'(x)fi - 4b'(x)V(x) = 0, (69) 

which can remarkably be converted to an EP equation, upon the scaling b(x) = 6 1 / 2 (x). 
Then, combining the knowledge of solvable EP cases (as per the discussion above) 
with that of the spatial profiles of the various (plane wave, solitary wave and elliptic 
function) solutions of the standard NLS equation for U, we can obtain special cases 
of g(x) for which explicit analytical solutions are available [211]. 



4-2. Methods from the linear limit 

When considering the GP equation as a perturbation problem, one way to do so 
is to consider the underlying linear Schrodinger problem, obtain its eigenvalues and 
eigenfunctions; subsequently one can consider the cubic nonlinear term within the 
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realm of Lyapunov-Schmidt (LS) theory (see Chap. 7 in Ref. [212]), to identify the 
nonlinear solutions bifurcating from the linear limit. 

We will discuss this approach in a general ID problem, with both a magnetic trap 
and an optical lattice potential, 

V(x) = Vmt(x) + Vol(x) = ifiV + V cos(2x), (70) 

following the approach of Ref. [213]. Considering the linear problem of the GP 
equation, using ip(x, t) — exp(—iEt)u(x) (where E is the linear eigenvalue) and 
rescaling spatial variables by CI 1 / 2 one obtains 

1 d 2 u 1 Vo / x \ E 

£t( = -n33 + ^ U+—COS 2-— )u=-u. (71) 



2 dx 2 2 Cl \ VI 1 / 2 ) Cl 

If we work in the physically relevant regime of < CI <C 1 and for Vo/Cl = 0(1), then 
one can use \x — Cl 1 / 2 as a small parameter and develop methods of multiple scales 
and homogenization techniques [213] in order to obtain analytical predictions for the 
linear spectrum. In particular, the one-dimensional eigenvalue problem using as fast 
variable X — x/fi and setting A = E^/Cl becomes: 

d 2 



^ CMT - fi dxdX +C ° L 



/Li Am, (72) 



where 



Id 2 1 

£MT= ~2W + 2 Xt (73) 

Col = -\^ + V cos(2X). (74) 

One can then use a formal series expansion (in fj,) for u and A 

u = u + y,u\ + n 2 u 2 + . . . , (75) 

A= A^ + A_i +Ao+ (76) 

Substitution of this expansion in the eigenvalue problem of Eq. (72), and after tedious 
but straightforward algebraic manipulations and use of solvability conditions for the 
first three orders of the expansion (0(1), 0(m) and 0(^ 2 )), yields the following result 
for the eigenvalue and the corresponding eigcnfunction of the eigenvalue problem of 
the original operator. The relevant eigenvalue of the n-th mode is approximated by: 

E^=-\v 2 + (l-\v 2 )n(n+ 1 ^, (77) 



and the corresponding eigenfunction is given by: 



u W (i) = c n H n ( X e ^¥ x ^ 



1 Vo ( 2x \ 

y cos yet}/ 2 ) 



(78) 



where c„ = (2"n!v / 7r) _ ^ 1//2 ^ is the normalization factor and H n (x) = 
e~ x (—l) n {d n /dx n )e x are the Hermite polynomials. 

Considering now the nonlinear problem Cu = —su 3 through LS theory, we obtain 
the bifurcation function 

G(p, AE) = -AEn + s ((u (n) ) 2 , (u (ll) ) 2 \ m 3 , (79) 
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Figure 1. (Color online) The first three (left to right) eigenfunctions of the 
linear Schrodinger equation with the potential of Eq. (70). The thick gray solid 
line corresponds to the eigenfunction for the case of VmtMi while the thin red 
solid and blue dashed lines correspond to the eigenfunction for Vmt(^) + Voh( x )t 
as computed numerically and theoretically [see Eq. (78)], respectively. The full 
potential VmtW + Vol(cg) (shifted for visibility, see scale on the right axis) is 
illustrated by the green dash-doted line. The parameters used in this example 
correspond to Vq =0.5 and Q = 0.1. 



for bifurcating solutions U n — fxu^ n ' , which bifurcate from the linear limit of E = E^ , 
with AE = E - E( n K The notation (/, g) = J f(x)g(x)dx will be used to denote the 
inner product of / and g. This calculation shows that a nontrivial solution exists 
only if sAE > (i.e., the branches bend to the left for attractive nonlinearities with 
s = — 1 and to the right for repulsive ones with s = 1) and the nonlinear solutions 
are created via a pitchfork bifurcation (given the symmetry u — > —u) from the linear 
solution. The bifurcation is subcritical for s = — 1 and supercritical for s = 1. Typical 
examples of the relevant solutions of the linear problem (from which the nonlinear 
states bifurcate) are shown in Fig. 1. 

Notice that for Vq — the problem becomes the linear quantum harmonic 
oscillator (parabolic potential) whose eigenvalues and eigenfunctions are known 
explicitly and are a special case of those of Eqs. (77)-(78) for Vq = 0. This perspective 
has been used in numerous studies as a starting point for numerical computation, e.g., 
in ID [105,214,215], as well as in higher dimensions [216]. 

It is natural to subsequently examine the stability of the ensuing nonlinear states, 
stemming from the linear problem. This can be done through the linearization of the 
problem around the nonlinear continuation of the solutions «W, with a perturbation 
u = w + iv. Then, the ensuing linearized equations for the eigenvalue A and the 
eigenfunction u can be written in the standard (for NLS equations) L + , L- form: 



L + w = 



L-v 



C + 3s (u^A 



w = -Xv (80) 



v = Xu. (81) 



Then, define n(L) and z(L) the number of negative and zero eigenvalues respectively 
of operator L, k r , k~ and k c the number of eigenvalues with, respectively, real positive, 
imaginary with positive imaginary part and negative Krein sign (see below) and 
complex with positive real and imaginary part. The Krein signature of an eigenvalue 
A is the sign((w, L + w)). One can then use the recently proven theorem for general 
Hamiltonian systems of the NLS type of Ref. [217] based on the earlier work of 
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Refs. [218-221] (see also Ref. [222]) according to which: 

k r + 2fcr + 2k c = n(L + ) + n{LJ) - n(D), (82) 

where D = dN/dE (N is the number of atoms of the state of interest). Then 
from Sturm-Liouville theory and given that is the only eigenfunction of L- with 
eigenvalue (by construction), we obtain n(L-) — j and z(L-) = 1 for the eigenstate 
U = (j,u^ (since it possesses j nodes). Similarly, using the nature of the bifurcation, 
one obtains n(L + ) = j and n(D) = if s = 1, while n(L + ) = j + 1 and n(LJ) = 1 
if s = —1. Combining these results one has that k r + 2k^ + 2k c = 2j. More detailed 
considerations in the vicinity of the linear limit [213] in fact show that the resulting 
eigenvalues have to be simple and purely imaginary i.e., k r — k c = and, hence, 
each of the waves bifurcating from the linear limit is spectrally stable close to that 
limit. However, the nonlinear wave bifurcating from the j-th linear eigenstate has 
j eigenvalues with negative Krein signature which may result in instability if these 
collide with other eigenvalues (of opposite sign). That is to say, the state has j 
potentially unstable eigendirections. 

While the above results give a detailed handle on the stability of the structures 
near the linear limit, it is important to also quantify the bifurcations that may 
occur (which may also, in turn, alter the stability of the nonlinear states), as well 
as to examine the dynamics of the relevant waves further away from that limit. 
A reduction approach that may be used to address both of these issues is that of 
projecting the dynamics to a full basis of eigenmodes of the underlying linear operator. 
Notice that we have seen this method before in the reduction of the GP equation 
in the presence of a strong OL to the DNLS equation. Considering the problem 
iut — Cu + s|u| 2 u — lou, we can use the decomposition u{x) — ^2^LoCj(t)qj(x), 
where qj(x) arc the orthonormalized cigenstates of the linear operator C. Setting 
a kim = (lkQiQm,qj) and following Ref. [223] straightforwardly yields 

M 

iCj = (p 3 -Lj)c 3 +S a llm C kClC* m , (83) 

k,l,m— 

where are the corresponding eigenvalues of the eigenstates Qj(x). It is interesting 
to note that this system with M — > oo is equivalent to the original dynamical 
system, but is practically considered for finite M, constituting a Galerkin truncation 
of the original GP PDE. This system preserves both the Hamiltonian structure of 
the original equation, as well as additional conservation laws such as the L 2 norm 

IMIi. = ££oM a - 

A relevant question is then how many modes one should consider to obtain a 
useful/interesting/faithful description of the original infinite dimensional dynamical 
system. The answer, naturally, depends on the form of the potential. The above 
reduction has been extremely successful in tackling double well potentials in BECs 
[223-225], as well as in optical systems [226]. In this simplest case, a two-mode 
description is sufficient to extract the prototypical dynamics of the system with M = 2. 
Then the relevant dynamical equations become: 

ico = (Mo - ^)c + sa° 0Q \c \ 2 c + sa? 10 (2|c 1 | 2 c + c\c* Q ), (84) 
ici = (pi - w)ci + sa in |ci| 2 ci + sa? 10 (2|c | 2 ci + clc{), (85) 

where we have assumed a symmetric double well potential so that terms such as aj 00 or 
a" n disappear and a\ ol — a\ w = a\ w = . . . = J q^q\dx. None of these assumptions 
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is binding and the general cases have in fact been treated in Refs. [223,224]. An 
angle-action variable decomposition Cj = Pje 1 ^ 3 leads to 

Po = sa?ioPiPosin(2<ty), (86) 

8(f) = - Ap + s(a° OQO pl - a\ llP \) 

-sa° U0 (2 + 2co S (S(f)))(p 2 -pi), (87) 

where 8(f) is the relative phase between the modes. Straightforwardly analyzing the 
ensuing equations [particularly Eq. (86)], we observe that the nonlinear problem can 
support states with pi = and po ^ (symmetric ones, respecting the symmetry 
of the ground state of the double well potential). It can also support ones with 
po = and p\ ^ (antisymmetric ones); these states are not a surprise since they 
did exist even at the linear limit. However, in addition to these two, the nonlinear 
problem can support states with p ^ and P \ ^ 0, provided that sin(2<50) = 0. 
These are asymmetric states that have to bifurcate because of the presence of 
nonlinearity. A more detailed study of the second equation shows that, typically, 
the bifurcation occurs for ||u||| 2 > N c — A^/(3a° 10 — ao 00 ) for s = — 1 (attractive 
case) and is a bifurcation from the symmetric ground state branch, while it happens 
for ||w||| 2 > N c = A/x/(3a5 10 — a{ n ) in the s = 1 (repulsive case) and is a bifurcation 
from the antisymmetric first excited state (see also Fig. 2 which shows the relevant 
states and bifurcation diagram). Notice that this is a pitchfork bifurcation (since there 
are two asymmetric states born at the critical point, each having principally support 
over each of the two wells). It should be mentioned that although this bifurcation 
is established for the two-mode reduction, it has been systematically confirmed by 
numerical analysis of the GP PDE in Refs. [223,224] for the case of a magnetic trap 
and an optical lattice or a magnetic trap and a defect respectively and it has been 
rigorously proved for a decaying at infinity double well potential in Ref. [225] (in 
the latter the corrections to the above mentioned N c were estimated). Based on the 
nature of the bifurcation (but this can also be proved within the two-mode reduction 
and rigorously from the GP equation) , we expect the ensuing asymmetric solutions to 
be stable, destabilizing the branch from which they are stemming, as is confirmed in 
the numerical linear stability results of Fig. 2. It is also worthwhile to point out that 
such predictions (e.g., the stabilization of an asymmetric state beyond a critical power) 
have been directly confirmed in optical experiments [226] in photorefractive crystals, 
and also have a direct bearing on relevant BEC experiments analyzed in Ref. [75]. 
We note in passing that in the physics literature, the problem is often tackled using 
wavefunctions that are localized in each of the wells of the double well potential (as 
linear combinations of the states q\ and q 2 used herein) [227-229], especially to study 
Josephson tunneling (but also to examine self-trapping) [75] . We refer the interested 
reader to these publications for more details. 

Such a few-mode approximation has also been successfully used in the case of 
three wells (in connection to applications in experiments in photorefractive crystals) in 
Ref. [230] . Naturally in that case, one uses three modes for the relevant decomposition 
and the corresponding dynamical equations grow in complexity very rapidly (as 
numerous additional overlap terms become relevant and the analysis becomes almost 
intractable). Similarly, such Galerkin approaches can also be used in the case where 
there is only a magnetic trap, in which case the underlying basis of expansion 
becomes that of the Hermite-Gauss polynomials [231]. In that case, in addition to the 
persistence of the linear states and a detailed quantitative analysis of their linearization 
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Figure 2. (Color online) The left panels illustrate a typical double well problem 
in the focusing case, while the right ones in the defocusing case. The top row 
shows the squared L 2 norm of the solutions (TV) as a function of the eigenvalue 
parameter /i (illustrating in each case the bifurcation of a new asymmetric branch). 
The second row shows the instability of the solid (blue) branch (symmetric in the 
left and antisymmetric in the right) past the critical point through the appearance 
of a real eigenvalue. The third row shows a particular example of the profiles for 
each branch and the fourth row shows the spectral plane of the linearization 
around them (absence of a real eigenvalue indicates stability). 



spectrum that becomes available near the linear limit, one can importantly predict 
the formation of new types of solutions. An example of this type consists of the space- 
localized, time-periodic (i.e., breathing) solutions in the neighborhood of, e.g., the 
first excited state (which has the form of a dark soliton) of the harmonically confined 
linear problem [231]. 

Finally, we indicate that such methods from the linear limit can equally 
straightforwardly be applied to higher dimensional settings and be used to extract 
complex nonlinear states. For instance, considering the problem 

d 2 u 1 du 1 d 2 u 9 , ,9 , , 

where also a rotational stirring term (frequent to condensate experiments [52,53]) is 
included, one can use a decomposition into the linear states q m ,i(r) exp(U9) [232], 
where m is the number of nodes of q m ,i- Then the underlying linear problem has 
eigenvalues X m .i = 2(\l\ + 1) + 4m + Itt and e.g., for solutions bifurcating from A = 6, 
one can write 

u = (xiqx t o{r) + yxqo,i< cos(Z'0) + iy 2 qo,i< sin(Z'0)) e 1/2 , (89) 
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Figure 3. (Color online) A typical ring solution (top left), multipole solution 
(top right), soliton necklace (bottom left) and vortex necklace (bottom right) 
that can be obtained from the near-linear analysis of the 2D problem through 
Eqs. (90)— (92). Reprinted from Ref. [232] with permission. 



together with lu = 6 + e5u>, and derive algebraic equations for xi, y\ and y 2 '- 

= Xl [n + 2x\ + 9l {2\ yi \ 2 + yi + yl)] , 

3 1 

- c g m+92xt(2yi + yt) + -\yi\ 2 yi + |J/|(2yi - y{), 



vl) 



3 2 
4 y2 



(90) 
(91) 

(92) 



where /i = -sSuj / (g Q ir), g x = g Q .i'/go, 92 = go,i>/gi> and c g = g 2 /gi and g x = / rqf fi dr, 
gi> = J rq^ ydr, go.c — / r Qi oQq udr. From these equations one can find real solutions 
containing only x\ (ring solutions), only y\ (multipole solutions), both £i and y\ 
(soliton necklaces), as well as complex solutions also involving y 2 ^ such as vortices 
and vortex necklaces. A sampler of these solutions is illustrated in Fig. 3; more details 
can be found in Ref. [232] , where the stability of such states is also analyzed, leading to 
the conclusion that the most robust among them are the (soliton and vortex) necklace 
and the vortex states. 

4-3. Methods from the nonlinear limit 

We partition our consideration of such methods to ones that tackle the stationary 
problem (in connection to the existence and the stability of the solutions) and to ones 
that address the dynamics of the perturbed solitary waves. 



If.. 3.1. Existence and stability methods. Consider a general Hamiltonian system of 
the form: 

§ = JE'(v), (93) 

where the J matrix has the standard symplectic structure (J 2 = —I) and E = 
j {l/2)\\v T \ 2 + s\v\ 4 ]dx for the case of the GP equation without external potential 
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(although the formalism presented below is very general [217,233,234]). Given that 
the above model in the case of the GP equation has certain invariances (e.g., 
with respect to translation and phase shift), one can use the generator T^ of the 
corresponding semigroup T(exp(wt)) of the relevant symmetry to make a change of 
variables v(t) = T(exp(u>t))u(t), which in turn leads to du/dt = JEq(u;u>), where 
E' (u;ui) — E'(u) — J^ 1 T t ^u. Defining then the appropriate conserved functional 
Qui = (1/2) (J~ 1 T LLl u, u), we note that relative equilibria satisfying E' (u;w) — will 
be critical points of Eq(u; oj) — E(u) — Q u (u). Then the linearization problem around 
such a stationary wave Uq reads: 

JE%(u ;uj)w = Xw. (94) 

Given the symmetries of the problem, this linearization operator has a non-vanishing 
kernel since: 

JE' ( ;(u ;w)T U!i u o = 0, (95) 

JE%(uo;w)d Ui uo=T Ui uo, (96) 

where each i corresponds to one of the relevant symmetries and the latter equation 
provides the generalized eigenvectors of the operator. 

The consideration of the perturbed Hamiltonian problem with a Hamiltonian 
perturbation such that the perturbed energy is Eq(u) + eE\(u) was considered in 
Refs. [217,222,233,234] and a number of conclusions were reached regarding the 
existence and stability of the ensuing solitary waves. Firstly, a necessary condition for 
the persistence of the wave is: 

(E[(u ;uj),T Ui u ) =0, (97) 

for all i pertaining to the original symmetries. This is a rather natural condition 
intuitively since it implies that the perturbed wave is a stationary solution if it is a 
critical point of the perturbation energy functional. The condition is also sufficient 
if the number of zero eigenvalues z(M) of the matrix M t j — (T Wi u , E"(uq; w)T Wj w ) 
is given byn-fc s , where n is the multiplicity of the original symmetries and k s the 
number of symmetries broken by the perturbation. 

As a result of the perturbation, 2k s eigenvalues (corresponding to the k s broken 
symmetries) will leave the origin, and can be tracked by the following result proved 
by means of LS reductions in Rcf. [217]. The eigenvalues will be A = ^/e\\ + 0(e), 
where the correction Ai is given by the generalized eigenvalue problem: 

(D A? +M 1 )v = 0, (98) 

where the matrix of symmetries (Dq)^ = {d Ui UQ, Eq (uq; w)d Uj UQ). 

In addition to this perturbative result on the eigenvalues, one can obtain a general 
count on the number of unstable eigendircctions of a Hamiltonian system [217], using 
the functional analytic framework of Refs. [218-221] (see also Ref. [222] for a different 
approach). In particular, for a linearization operator C u — E"(u) — J^T^ and a 
symmetry matrix Dij — (<9 Wi u, 

k r + 2fcr + 2k c = n{C u ) - n(D) - z(D), (99) 

where the relevant symbolism has been introduced in Sec. 4.2. In fact, the latter 
subsection constitutes a special case example of this formula, in the case of the form 
of 
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We now give a special case example of the application of the theory in the presence 
of a linear and a nonlinear lattice potential of the form [235]: 

iu t = ~u xx - (1 + en 1 (x)) \u\ 2 u + en 2 (x)u. (101) 
Then the problem can be rephrased in the above formalism with 

Ei(u) = J + (n 2 (x)\u\ 2 - ^ ni {x)\u\ 4 ^ dx. (102) 

Therefore, as indicated above, the persistence of the stationary bright solitary wave of 
the form u = y/Jl scch[ v ^Z(x — S,)]e lS (with S = fi/2) is tantamount to: V^Ei(u) = 0. 
This implies that the wave is going to persist only if centered at the parameter-selected 
extrema of the energy (which are now going to form, at best, a countably infinite set 
of solutions, as opposed to the one-parameter infinity of solutions previously allowed 
by the translational invariance). 

Furthermore, the stability of the perturbed wave is determined by the location of 
the eigenvalues associated with the translational invariance; previously, the relevant 
eigenvalue pair was located at the origin A = of the spectral plane of eigenvalues 
A = A r + On the other hand, we expect the eigenvalues associated with the U(l) 
invariance (i.e., the phase invariance associated with the L 2 conservation) to remain 
at the origin, given the preservation of the latter symmetry under the perturbations 
considered herein. Adapting the framework of Ref. [217], we have that the matrices 
that arise in Eq. (98) are given by: 

f (d x u ,-xu ) o A / /iV2 o A 

D °-{ 2(u ,d,u ) )-{ o —yr x l 2 J' (103) 

and 

= ^/(^(«°) 2 -i^M 4 )^ oj (104) 

One can then use the above along with Eq. (98) to obtain the relevant translational 
eigenvalue as: 



/iV2 J \2 dx 2 y VJ 4 dx 2 

Based on this expression, the corresponding eigenvalue can be directly evaluated, 
provided that the extrema of the effective energy landscape E\ are evaluated first. 
This effective energy landscape V e s (£) = eE\ is a function of the solitary wave location 
£. The physical intuition of the above results is that the stability or instability of the 
configuration will be associated with the convexity or concavity, respectively, of this 
effective energy landscape. Some examples of the accuracy of such a prediction arc 
provided in Fig. 4, for specific forms of n\(x) and 712(2). 

This class of techniques has been applied to different problems with spatial 
variation of the linear [236] or nonlinear [237] potential. They can also be applied to 
multi-component problems [233] or to problems with different nonlinearity exponents 
[238] or higher dimensions [239]. We note in passing that in addition to these methods, 
for periodic variations of the potential, and for appropriate regimes (for details see 
Refs. [238-240]), one can develop multiple-scale techniques exploiting the disparity in 
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Figure 4. (Color online) Typical examples of the translational eigenvalue 
as obtained numerically (solid/blue line) versus the analytical prediction 
(dashed/green line). The linear and nonlinear potentials are: ni(x) = Acos(kix) 
and ri2(x) = Bcos(k2X + A</>). In the left panels we assume that A = B = 1, and 
fix &2 = 2-7r/5 and A<p = and examine the relevant translational eigenvalue (its 
real part and its square) as a function of k\. Notice that there is a transition from 
instability to stability as fci is increased. In the right panels we select A = B = 1 
and fci = fc2 = 27r/5 and vary A0 £ [0, 27r], Notice that in the latter case the 
(critical point associated with the stationary) location of the solitary wave also 
changes with A<f> and its theoretical and numerical values are also given (again in 
dashed and solid lines, respectively). Notice in all the cases the accuracy of the 
theoretical prediction. 



spatial scales between the solution and the potential. We refer the interested reader to 
the above references for further technical details. This type of averaging techniques is 
popular not only when the linear or nonlinear potential presents spatial variations of 
a characteristic scale, but also similarly when these variations are temporal [241-243], 
especially because it is more straightforward in the averaged equations to extract 
conclusions about the possible existence or potential collapse or dispersion of the 
solutions [238,239,244]. 

A similar approach can be used in the case of dark solitons in examining the 
persistence and stability of the waves in the presence of external potentials; however 
in the latter case, it is a much harder task to control the linearization spectrum of the 
problem since it encompasses the origin. This complication has allowed this problem 
to be tackled only recently at the nonlinear limit [245] and perturbatively away from 
that limit [246] . The main result of Ref. [245] is that by using the limit 

lira ,g(A) = lim ((L_A) _1 u', tt'), (106) 

one can infer the stability of the black soliton, since if this quantity is positive the 
soliton will have a real eigenvalue and will be unstable, while if non-positive, it will be 
stable. However, one of the problems with this expression is that even when the soliton 
is analytically available, it is relatively hard to evaluate (see, e.g., the example of the 
integrable cubic case worked out in Ref. [245]). On the other hand, although there 
exist results on the orbital stability of dark solitons [247] and of other structures such 
as bubbles (black solitons with zero phase shift) [248] (see Ref. [246] for a more detailed 
discussion of earlier works), the work of Ref. [246] was the first one to establish detailed 
estimates on the relevant eigenvalues, using once again the technique of Lyapunov- 
Schmidt reductions in the limit of small potential perturbations. The main results 
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can be summarized as follows. For the equation: 
iu t = ~u xx + f(\u\ 2 )u + eV{x)u 
(i) The analogous condition to the persistence condition (97) is now: 
M'(0 - / V'(x)[q - ul(x - 0}dx = 0, 



(107) 



(108) 



where the unperturbed solution uq asymptotes to i^/qa at ±oo; i.e., the 
background of the solution is now appropriately incorporated in Eq. (108) [in 
comparison to Eq. (97)] . 

(ii) The dark soliton will be spectrally unstable in the GP equation with exactly one 
real eigenvalue (for small e) in the case where 



while it will be unstable due to two complex-conjugate eigenvalues with positive 
real part if M"(£) > 0. This is the analogous condition to the curvature 
of the effective potential; however, notice the disparity of this condition from 
what would be expected intuitively based on the notion of convexity /concavity. 
The latter result is a direct byproduct of the nature of the essential spectrum 
(encompassing the origin in this defocusing case), which upon bifurcation of the 
translational eigenvalue along the imaginary axis makes it directly complex (case 
of M"(£) > 0). The location of the relevant eigenvalue in the GP case of cubic 
nonlincarity is given to leading order by: 



which is consonant with the above result. A form of this expression for general 
nonlincarities was also obtained in Rcf. [246] . 

A case example of the possible dark soliton solutions for a potential of the form 
V(x) = x 2 exp(— k\x\) is shown in Fig. 5, illustrating the quantitative accuracy 
of Eqs. (108)-(110). In this case, the formalism elucidates a subcritical pitchfork 
bifurcation whereby three dark soliton solutions (one centered at £ = and two 
symmetrically at £ ^ 0) eventually merge into a single unstable kink centered at 



One of the fundamental limitations of this result is its being dependent upon the 
decaying nature of the potential at ±oo. The fundamentally different nature of the 
spectrum in the presence of parabolic or periodic potentials (or both) makes it much 
harder to provide such considerations in the latter cases. While this can be done in 
some special limits (such as the Thomas-Fermi, large chemical potential limit of the 
appendix of Ref. [231]), in that setting it is generally easier to use the linear limit 
approach of Sec. 4.2. 

4-3.2. Perturbation theory for solitons. Dynamics of either bright or dark matter- 
wave solitons can be studied by means of the perturbation theory for solitons [249-252] 
(see also Ref. [253] for a review). Here, we will briefly discuss an application of this 
approach upon considering the example of soliton dynamics in BECs confined in an 
external potential, say V(x), which is smooth and slowly-varying on the soliton scale. 
This means that in the case, e.g., of the conventional parabolic trap V(x) = (l/2)£l 2 x 2 , 




(109) 




(110) 
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Figure 5. (Color online) An example of a subcritical pitchfork bifurcation in the 
parameter k of the potential V{x) = x 2 exp(— k\x\) for fixed e = 0.2 in Eq. (107). 
The left panel shows the center of mass so = £ of the dark soliton kink modes 
(so ^ by dashed line, so = by thick solid and dashed lines). The theoretical 
predictions of so based on Eq. (108) are shown by dash-dotted line. The vertical 
line gives the theoretical prediction for the bifurcation point k = reo = 3.21. The 
right panel shows the real part of the unstable eigenvalues for the corresponding 
solutions, using the same symbolism as the left panel. The theoretical predictions 
of eigenvalues are shown by thick and thin dash-dotted lines, respectively for the 
branches with so = and so ^ 0. Notice that for the quartet of eigenvalues of 
the branch centered at the origin, the small jumps are due to the finite size of the 
computational domain (see Ref. [246] for details). 




the effective trap strength is taken to be fl ~ e, where £ < 1 is a formal small 
(perturbation) parameter. Taking into regard the above, we consider the following 
perturbed NLS equation 

du 1 d 2 u . . . 

with the perturbation being the potential term R(u) = V(x)u, and g = ±1 
corresponding to repulsive and attractive interactions. Soliton dynamics in the 
framework of Eq. (Ill) can then be treated perturbatively, assuming that a perturbed 
soliton solution can be expressed as 

u(x, t) = u s (x, t) + eud(x, t) + eu r (x, t). (112) 

Here, u s (x,t) has the same functional form as the soliton solutions (24) for g = — 1 
and (27) for g = +1, but with the soliton parameters depending on time. On the 
other hand, Ud is a function localized near the soliton, describing the deformation 
of the soliton (i.e., the change of the soliton shape) and u r is the radiation (in the 
form of sound waves) emitted by the soliton. In fact, the effect described by Ud is 
not significant, as the small change in the soliton shape does not modify its motion, 
while the emission of radiation may be neglected for sufficiently weak perturbations. 
Thus, here we will consider solely the first term in Eq. (112), which corresponds to 
the so-called adiabatic approximation of the perturbation theory for solitons [253]. 

First we discuss the dynamics of bright solitons in external potentials. Taking 
into regard that for g = —1 and R(u) = 0, Eq. (Ill) has a bright soliton solution of 
the form given in Eq. (24), we assume that in the perturbed case with R(u) ^ a 
soliton solution can be expressed as 

u(x, t) = t\ sech[r;(x — 2:0)] enp[i(kx — <fi(t)] (H3) 

where xq is the soliton center, the parameter k = dxo/dt defines both the soliton 
wavenumber and velocity, and 4>(t) = (l/2)(fc 2 — i] 2 )t is the soliton phase. In the 
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case under consideration, the soliton parameters 77, k and x are considered to be 
unknown, slowly-varying functions of time t. Then, from Eq. (Ill), it is found that 
the number of atoms N and the momentum P (which are integrals of motion of the 
unperturbed system), evolve, in the presence of the perturbation, according to the 
following equations, 



= -21m 

dt 



r+00 ip r r-+oc o * 

/ R{u)u*dx , — = 2Rc / R(u)^-dx 
J-00 at [J_ oc dx 



(114) 



We now substitute the ansatz (113) (but with the soliton parameters being functions 
of time) into Eqs. (114) and obtain the evolution equations for rj(t) and k(t), 

dn „ dk dU 

i = °> dt=~^ ^ 

where U(x n ) is given by the expression of Eq. (55). In the case, however, of slow 
variation of the potential on the scale of the solitary wave (the case of interest here) , 
a simple Taylor expansion yields the same equation but with U = V i.e., the trapping 
potential. 

Now, recalling that dxo/dt = k, we may combine the above Eqs. (115) to derive 
the following equation of motion for the soliton center: 

d 2 *o dV 

(116) 

which shows that the bright matter-wave soliton behaves effectively like a Newtonian 
particle. Note that in the case of a parabolic trapping potential, i.e., V(x) = 
(l/2)Q 2 x 2 , Eq. (116) implies that the frequency of oscillation is Q; this result is 
consistent with the Ehrcnfcst theorem of the quantum mechanics, or the so-called 
Kohn theorem [198], implying that the motion of the center of mass of a cloud of 
particles trapped in a parabolic potential is decoupled from the internal excitations. 
Note that the result of Eq. (116) can be obtained by other methods, such as the WKB 
approximation [254], or other perturbative techniques [255-257]. 

We now turn to the dynamics of dark matter-wave solitons in the framework 
of Eq. (Ill) for g = +1. First, the background wavefunction is sought in the form 
u = $(x) exp(— i[it) (/1 being the normalized chemical potential), where the unknown 
function &(x) satisfies the following real equation, 

^ + l^-& = V(x)$. (117) 

Then, following the analysis of Ref. [258], we seek for a dark soliton solution of 
Eq. (Ill) on top of the inhomogeneous background satisfying Eq. (117), namely, 
u = $(#) exp(— i/j,t)v(x, t), where the unknown wavefunction v(x, t) represents a dark 
soliton. This way, employing Eq. (117), the following evolution equation for the dark 
soliton wave function is readily obtained: 

»7* + 2ft?-*>' - 1)v = -T x l ^dx- (118) 

Taking into regard the fact that in the framework of the Thomas- Fermi approximation 
the profile can be simply approximated by Eq. (117), Eq. (118) can be simplified to 
the following dcfocusing perturbed NLS equation, 

l % + Y^~ M2 - l)v = Q{v) ' (119) 
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with the perturbation Q(v) being of the form, 



,. 1 dV dv 



= (1 - «V + 5( _ T ^__. (120) 

In the absence of the perturbation, Eq. (119) has the dark soliton solution (for fi = 1) 
v(x,t) = cos<ptanh£ + ismp 1 where Q = cos [x — (sin ip)t] (recall that cos<p and 
sin if are the soliton amplitude and velocity respectively and ip is the soliton phase 
angle). To treat analytically the effect of the perturbation (120) on the dark soliton, 
we employ the adiabatic perturbation theory devised in Ref. [114]. Assuming, as in 
the case of bright solitons, that the dark soliton parameters become slowly-varying 
unknown functions of t, the soliton phase angle becomes ip — ► ip(t) and, as a result, the 
soliton coordinate becomes C - * C = cosip(t) [x — Xo(t)], where Xo(t) = j Q sin (p(t')dt' , 
is the soliton center. Then, the evolution of the parameter ip is governed by [114], 

dp 1 



-Re 



dv* ' 
Q(v)—dx 



(121) 



dt 2 cos 2 ip sin tp 
which, in turn, yields for the perturbation in Eq. (120): 
d<p 1 dV 

i = -2 C °^8x- - (122) 

To this end, combining Eq. (122) with the definition of the dark soliton center, we 
obtain the following equation of motion (for nearly stationary dark solitons with 
cos (p w 1), 

dt* - 2dx - [iZ6) 

Equation (123) implies that the dark soliton, similarly to the bright one, behaves 
like a Newtonian particle. However, in an harmonic potential with strength f2, the 
dark soliton oscillates with another frequency, namely f2/\/2, a result that may be 
considered as the Ehrenfest theorem for dark solitons. The oscillations of dark solitons 
in trapped BECs has been a subject that has attracted much interest; in fact, many 
relevant analytical works have been devoted to this subject, in which various different 
perturbative approaches have been employed [259-263]. It should also be mentioned 
that a more general Newtonian equation of motion, similar to Eq. (123) but also valid 
for a wider class of confining potentials, was recently discussed in Ref. [246] . 

Perturbation theory for dark solitons may also be applied for dark solitons with 
radial symmetry, i.e., for ring or spherical dark solitons described by a GP equation 
of the form 

;^ = -^vV + IVO + vw, (124) 

where V 2 = <9 2 + (D — l)r~ 1 d r is the transverse Laplacian, V(r) = (l/2)fi 2 r 2 , and 
D = 2, 3 correspond to the cylindrical and spherical case, respectively. In this case, 
Eq. (124) can be treated as a perturbed ID defocusing NLS equation provided that 
the potential term and the term <~ r _1 can be considered as small perturbations; 
this case is physically relevant for weak trapping potentials with <C 1, and radially 
symmetric solitons of large radius r . Then, it can be found [264] (see also Ref. [262]) 
that the radius r of the radially symmetric dark solitons is governed by the following 
Newtonian equation of motion, 

IF " -<W- < 125 > 
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where the effective potential is given by V e g(ro) = (l/4)£l 2 r 2 — lnr^ 15 [note 
that in the ID limit of D = 1, Eq. (125) is reduced to Eq. (123)]. Clearly, in 
this higher-dimensional setting the equation of the soliton motion becomes nonlinear, 
even for nearly black solitons, due to the presence of the repulsive curvature-induced 
logarithmic potential. We finally note that such radially symmetric solitons are 
generally found to be unstable, as they either decay to radiation (the small-amplitude 
ones) or are subject to the snaking instability (the moderate- and large-amplitude 
ones) , giving rise to the formation of vortex necklaces [264] . 



4-3.3. The reductive perturbation method. Another useful tool in the analysis of 
the dynamics of matter-wave solitons (and especially the dark ones) is the so-called 
reductive perturbation method (RPM) [265]. Applying this asymptotic method, one 
usually introduces proper "stretched" (slow) variables to show that small- amplitude 
nonlinear structures governed by a specific nonlinear evolution equation can effectively 
be described by another equation. Such a formal connection between soliton equations 
was demonstrated for integrable systems in Ref. [266], and then extended to the 
reduction of nonintegrable models to integrable ones, first in applications in optics (see, 
e.g., Refs. [267-269]) and later in BECs. Here, we will briefly describe this method 
upon considering, as an example, an inhomogeneous generalized NLS equation similar 
to Eq. (107), namely, 

iu t = -^u xx + V(X)u + g(p)u, (126) 

which is characterized by a general nonlinearity g(p) (depending on the density 
p = \u\ 2 ), and a slowly varying external potential V(X), depending on a slow variable 
X = e 3 / 2 x (with e being a formal small parameter). Our main purpose is to show that 
this rather general NLS-type mean-field model can be reduced to the much simpler 
Korteweg-dc Vrics (KdV) equation with variable coefficients. The latter, has been 
used in the past to describe shallow water-waves over variable depth, or ion-acoustic 
solitons in inhomogeneous plasmas [270] , and, as we will discuss below, it provides an 
effective description of the dark soliton dynamics in BECs. 

Following the analysis of Ref. [150], we first derive from Eq. (126) hydrodynamic 
equations for the density p and the phase 0, arising from the Madelung transformation 
u = y/pex-p(i(j)), and then introduce the following asymptotic expansions, 

p = p (X) + e Pl (X,T) + e 2 P2 (X,T) + ■ ■ ■ , (127) 

cj> = - p t + e 1/2 ^i (X, T) + e 3 / 2 2 (X, T) + ■ ■ ■ , (128) 

where po(X) is the ground state of the system determined by the Thomas- Fermi 
approximation g(po) = Po — V(X) (with p being the chemical potential), and 
T = e 1 / 2 (t — J* C~ 1 (x')dx') is a slow time-variable [where C = \fgaPo is the local 
speed of sound and go = {dg / dp)\ p=Pa \. This way, in the lowest-order approximation 
in e we obtain an equation for the phase, 

MX,T) = -g (X) ( T Pl (X,T')dT', (129) 
Jo 



and derive the following KdV equation for the density pi, 



(3g + po'go) 1 d 

P1X 2C* PlPlT + 8C^ P1TTT ~ ~dX 



ln(|C|<7o) 1/2 l Pi, (130) 
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where g = (d 2 g/dp 2 )\ p=po . Importantly, in a homogeneous gas with pa(X) = p p = 
const., Eq. (130) is the completely integrable KdV equation; the soliton solution of the 
latter, which is ~ sech 2 (T), is in fact a density notch on the background density p p [see 
Eq. (127)], with a phase jump across it [see Eq. (129), which implies that (f>\ ~ tanh(T)] 
and, thus, it represents an approximate dark soliton solution of Eq. (126). On the other 
hand, the results obtained earlier for the analysis of the KdV equation with variable 
coefficients [271,272] have been used to analyze shallow soliton dynamics in the BEC 
context. In particular, the KdV Eq. (130) was obtained in the framework of the cubic 
nonlinear version of Eq. (126) (with g(p) <~ p), and used to study the dynamics and 
the collisions of shallow dark solitons in BECs [260,273]. Moreover, other versions of 
Eq. (130) relevant to the Tonks gas (corresponding to g(p) ~ p 2 ) [159], as well as the 
BEC-Tonks crossover regime (corresponding to a generalized nonlinearity) [150] were 
derived and analyzed as well. 

Finally, it is relevant to note that there exist studies in higher-dimensional (disk- 
shaped) BECs, where the RPM was used to predict 2D nonlinear structures, such 
as "lumps" described by an effective Kadomtsev-Petviashvili equation [274], and 
"dromions" described by an effective Davey-Steartson equation [275,276]. 

4-4- Methods for discrete systems 

As our final class of methods, we will present a series of results that are relevant 
to systems with periodic potentials and, in particular, with discrete lattices per the 
Wannier function reduction of Sec. 3.6.2. Starting with the prototypical discrete model 
of the ID DNLS equation (see, e.g., Ref. [186] for a review) of the form: 

iu n = -e(u n+ i + u„_i) - \u n \ 2 u n , (131) 

we look for standing waves of the form: u n = exp(ipt)v n which satisfy the steady 
state equation 

(M- \v n \ 2 )v n = e{v n +i +v n - 1 ). (132) 

One of the fundamental ideas that we exploit in this setting is the so-called anti- 
continuum (AC) limit of MacKay and Aubry [277] for e = 0, where Eq. (132) is 
completely solvable v n = {0, ± v /7iexp(i6' n )}, where 9 n is a free phase parameter for 
each site. However, a key question is which ones of all these possible sequences of 
v n will persist as solutions when e =/= 0. A simple way to see this in the ID case of 
Eq. (131) is to multiply Eq. (132) by w* and subtract the complex conjugate of the 
resulting equation, which in turn leads to: 

v n v n+i - v n Vn+i = const. 2arg(w„ + i) = 2arg(w„), (133) 

since we are considering solutions vanishing as in ±oo. Without loss of generality 
(using the scaling of the equation), we can scale p, = 1, in which case the only states 
that will persist for finite e are ones containing sequences with combinations of v n = ±1 
and v n = 0. A systematic computational classification of the simplest ones among 
these sequences and of their bifurcations is provided in Ref. [278]. Notice that we are 
tackling here the focusing case of s = — 1, however, the defocusing case of s — 1 can 
also be addressed based on the same considerations and using the so-called staggering 
transformation w n = (—l) n u n (which converts the defocusing nonlinearity into a 
focusing one, with an appropriate frequency rescaling which can be absorbed in a 
gauge transformation). 
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lim(6„,£; 1 6„) = -- =► A 2 = 2 7 = 2(b n ,£_6„) (137) 



We subsequently consider the issue of stability, using once again the standard 
symplectic formalism JCw = Au>, where C is given by Eq. (100) and J is the symplectic 
matrix. In this case, the L + and L_ operators are given by: 

(1 - 3v 2 )a„ - e(a„+i + o n _i) = £+a„ = -A6 n (134) 

(1 - O&n - e(&„+i + & n -i) =C-b n = \a n . (135) 
We again use the AC limit where we assume a sequence for w„ with A "excited" (i.e., 
7^ 0) sites; then, it is easy to see that for e = these sites correspond to eigenvalues 
\ L+ = —2 for L + and to eigenvalues X L _ = for L_, and they result in TV eigenvalue 
pairs with A 2 = for the full Hamiltonian problem. Hence, these eigenvalues are 
potential sources of instability, since for e ^ 0, A — 1 of those will become nonzero 
(there is only one symmetry, namely the U(l) invariance, persisting for e ^ 0). The 
key question for stability purposes is to identify the location of these N — 1 small 
eigenvalue pairs. One can manipulate Eqs. (134)-(135) into the form: 

C-b n = -A 2 £; 1 6„ =* A 2 = - (136) 

("71) •'-+ O n J 

Near the AC limit, the effect of L + is a multiplicative one (by —2). Hence: 

1 

;s^"'~+ ~ 7 " 2 

Therefore the problem reverts to the determination of the spectrum of L_. However, 
using the fact that v n is an eigenfunction of L_ with Xl_ = and the Sturm 
comparison theorem for difference operators [279], one infers that if the number 
of sign changes in the solution at the AC limit is m (i.e., the number of times 
that adjacent to a +1 is a —1 and next to a —1 is a +1), then n(LJ) = m 
and therefore from Eq. (137), the number of imaginary eigenvalue pairs of C is m, 
while that of real eigenvalue pairs is consequently (N — 1) — to. An immediate 
conclusion is that unless m — N — 1, or practically unless adjacent sites are out- 
of-phase with each other, the solution will be immediately unstable for e 0. 
Notice that this is also consistent with the eigenvalue count of Refs. [217,222] since 
n(C) - n(D) = (N + m) - 1 = (N + m - 1) + 2 x m + 2 x = k r + 2fcr + 2k c (it 
is straightforward to show by the definition of the Krein signature [280] to show that 
these m imaginary pairs have negative Krein signature). 

One can also use Eq. (137) quantitatively to identify the relevant eigenvalues 
perturbatively for the full problem by considering the perturbed (originally zero when 
e = 0) eigenvalues of L_ in the form: 

^n* 1 -^ (138) 

where £_ = + eC^ + 0{e 2 ) and a similar expansion has been used for the 
eigenvector b n . Also Xl_ = cyi + 0(e 2 ). Projecting the above equation to all 

the eigenvectors of zero eigenvalue of one can explicitly convert Eq. (138) 

into an eigenvalue problem [280] of the form Mc — jic, where the matrix M has 
off-diagonal entries: M n>n +i = M n+ i i7l = — cos(#„+i — 9 n ) and diagonal entries 
M B ,„ = (cos(# n _i — 6 n ) + cos(0„+i — n )). Then it is straightforward to compute 
71 and subsequently from Eq. (137) to evaluate the corresponding A = ±y/2eyi. For 
example, for one-dimensional configurations with two-adjacent sites with phases 0\ 
and 02, the matrix M becomes: 

M — ( ™<*"*> - C °ff d39) 
V - cos(0i - 6 2 ) cos(#i -62) J 



CONTENTS 



42 




Figure 6. (Color online) A typical example of a two-site configuration is shown in 
the left panels of the figure, and a corresponding one of a three-site configuration 
in the right panels. The top panels show a solution (for a particular value of e) 
and its corresponding spectral plane (A r , A;) of eigenvalues A = X r +i\i, while the 
bottom ones show the dependence of the relevant real eigenvalues as a function of 
the inter-site coupling e obtained analytically (dashed/red lines) and numerically 
(solid/black lines). 



whose straightforward calculation of eigenvalues leads to A 2 = and A 2 = 
±2-\/ecos(#i — 6*2). Notice that this result is consonant with our qualitative theory 
above since for same phase excitations (6± = 6 2 ), the configuration is unstable, 
while the opposite is true if B\ = 62 ± tt. Similar calculations are possible for 3- 
site configurations with phases #1,2,3, in which case, one of the eigenvalues of M is 
again (as has to generically be the case, due to the U(l) invariance), while the other 
two are given by: 

71 = cos(# 2 - 9i) + cos(6>3 - 9 2 ) 

± Vcos 2 (0 2 cos{e 2 - 9 X ) cos(# 3 - O2) + cos 2 (#3 - 2 ) . (140) 

Some of the examples of the accuracy of these theoretical predictions for some typical 
two-site and three-site configurations (in particular, the in-phase ones, which should 
have one and two real eigenvalue pairs respectively) are shown in Fig. 6. 

This approach can be generalized to different settings, such as in particular higher 
dimensions [281,282] or multi-component systems [283]. Perhaps the fundamental 
difference that arises in the higher dimensional settings is that one can excite sites 
over a contour and then, for the N excited sites around the contour the persistence 
(Lyapunov-Schmidt) conditions can be obtained as a generalization of Eq. (131) that 
reads [281,284]: 

sin(0i - 6 2 ) = sin(6? 2 -63) = ... = sm(9 N - 6 X ), (141) 

which indicates that a key difference of higher dimensional settings is that not 
only "solitary wave" structures with phases 8 E {0, tt} are possible, but also both 
symmetric and asymmetric vortex families [281,284] may, in principle, be possible 
[although Eq. (141) provides the leading order persistence condition and one would 
need to also verify the corresponding conditions to higher order to confirm that 
such solutions persist [281]]. Such vortex solutions had been predicted numerically 
earlier [285,286] and have been observed experimentally in the optical setting of 
photorefractive crystals [287,288]. Performing the stability analysis is possible for 
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these higher dimensional structures, although the relevant calculations are technically 
far more involved. However, the theory can be formulated in an entirely general 
manner: we give its outline and some prototypical examples of higher dimensional 
theory-computation comparisons below. 

The existence problem can be generally formulated in the multi-dimensional case 
as the vanishing of the vector field F„ of the form: 

(1- |0„| 2 )0„-e£0„ 
(1 - |0„| 2 )0* - e£</>* 
If we then define the matrix operator: 

1-2|0„| 2 -4>l 
-4>l i-2|<M 2 



(142) 



H n = 



— e { s +ei + s -ei + s +e 2 + s -e 2 + s +e 3 + s - 



1 
1 



(143) 



where the s± ei denotes the shift operators along the respective directions, then the 
stability problem is given by aTlij) = where the 2-block of a is the diagonal 

matrix of (1, —1) at each node n. Furthermore, the existence problem is connected to 
H through: H = D^F(<p, 0). At the AC limit of e = 



(W (0) )„ = 



1 
1 



, n£ S A 



(« (0) )n = 



-1 

-2i0„ 



, n e S, 



where S is the set of excited sites. Then the eigenvectors of zero eigenvalue will be of 
the form: 



(e„)k = i 



Defining the projection operator: 

(e„,f) 



( e >. 



'(fi)„-e ie "(f 2 )„) 



neS, (144) 
(145) 



and decomposing the solution as: 

= (o) (0) + v eX, 

one can obtain the Lyapunov-Schmidt persistence conditions as [282]: 

g(0, e) = 7>F(0 (O) (0) + <p(0, e), e) = 0. (146) 

This leads to the persistence theorem [282]: The configuration <f^ Q \9) can be 
continued to the domain e G O(0) if and only if there exists a root 0* of the vector 
field g(0, e). Moreover, if the root 0* is analytic in e e O(0) and 9* = 6 + O(e), the 
solution cf> of the difference equation is analytic in e £ O(0), such that 



(147) 



k=l 



One can also formulate on the same footing a general stability theory. More 
specifically, let the solution of interest persist for e ^ 0. If the operator W has a 
small eigenvalue /x of multiplicity d, such that /j = e k \i^ + 0(e fc+1 ), then the full 
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Hamiltonian eigenvalue problem admits (21?) small eigenvalues A. These are such 
that A = e k / 2 X k/2 + 0(e k / 2+1 ), where non-zero values A fe / 2 are found from 

odd k: M {k) a = lj\ 2 k/2 a, (148) 
event. M™a + \x k ,2^ k) oc = \\\ /2 a, (149) 

where 

£ (fe) = v[n^^ k ' ) (e ) + - + n^ +1 ^ (0 He ) 

and k' = (k — l)/2. For more details, we refer the interested reader to Ref. [282]. 

In Fig. 7 we show typical examples of 2D and 3D configurations that satisfy 
the persistence conditions formulated above. The former configuration is a vortex of 
topological charge 5 = 2 (i.e., its phase changes uniformly by n/2 from each site to the 
next so that it changes from to At: around the discrete contour of the solution) . This 
structure is unstable due to a real eigenvalue pair that is theoretically predicted from 
Eqs. (148)-(149) to be A = ±y/ \/80 — 8e (while it also has a pair of simple eigenvalues 

A = ±ieVv80 + 8, a quadruple eigenvalue A = ±iey/2 and a single eigenvalue of 
higher order). The latter configuration is a three-dimensional diamond configuration 
(a quadrupole in the plane with phases 0, n, 0, ir and two out-of-plane sites with phases 
7r/2 and 3tt/2). This is a stable 3D structure with a single eigenvalue A = ±4ie, a triple 
eigenvalue A = ±2ie and an eigenvalue of higher order according to Eqs. (148)-(149). 
Notice in both cases the remarkable agreement between the theoretical prediction for 
the eigenvalues (dashed lines) and the full numerical results (solid lines). 



5. Special topics of recent physical interest 

In this section we give a very brief overview of some of the more recent themes 
of interest in the nonlinear phenomenology emerging in the realm of Bosc-Einstcin 
condensation. We pay special attention to the physical motivation of the different 
topics. Note that an in-depth treatment of the emergent nonlinear behavior displayed 
by BECs and the synergy between experiments and theory can be found in the recent 
review [55]. 



5.1. Spinor/Multicomponent condensates 

Advances in trapping techniques for BECs have opened the possibility to 
simultaneously confine atomic clouds in different hyperfine spin states. The first 
such experiment, the so-called pseudospinor condensate, was achieved in mixtures of 
two magnetically trapped hyperfine states of 87 Rb [289]. Subsequently, experiments 
in optically trapped 23 Na [290] were able to produce multicomponent condensates 
for different Zccman sub-levels of the same hyperfine level, the so-called spinor 
condensates. In addition to these two classes of experiments, mixtures of two different 
species of condensates have been created by sympathetic cooling (i.e., condensing one 
species and allowing the other one to condense by taking advantage of the coupling 
with the first species) were 41 K atoms were sympathetically cooled by 87 Rb atoms 
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Figure 7. (Color online) The top panel shows a 2D S = 2 vortex configuration 
(left panels show its real and imaginary part -top- as well as amplitude and 
phase -bottom-) and its linear stability real and imaginary eigenvalues (right 
panels). Same thing in the bottom for a diamond 3D configuration (see the 
explanation for its phases in the text). The diamond configuration is shown 
using iso-level contours of different hues: blue/red (dark gray/gray in the black- 
and-white version) are positive/negative real iso-contours while the green/yellow 
(light gray/very light gray in the black-and-white version) correspond to 
positive/negative imaginary iso-contours. The theoretical predictions for the 
eigenvalues as a function of the coupling are shown by dashed line, while the 
full numerical results by solid ones. Partially reprinted from Ref. [280] with 
permission. 



[291]. More exotic mixtures are also being currently explored in degenerate fermion- 
boson mixtures in 40 K- 87 Rb [292] and pure degenerate fermion mixtures in 40 K- 6 Li 
[293-295]. 

The mean-field dynamics of such multicomponent condensates is described by a 
system of coupled GP equations analogous to Eq. (43), that for mixtures of AT bosonic 
components reads 

dtp n 1 ~ 

i-Qj- = --Alp n + V n (r)lp n + 2J [g n ,k\^k\ 2 ^n - K nt k4>k + A/Xfe^n] ~ i^n4>n, (150) 

were ip n is the wavefunction of the n-th component (n = 1,...,AT), V n (r) is the 
potential confining the n-th component, A/x^fc is the chemical potential difference 
between components n and k, and a n describes the losses of the n-th component. The 
components n and k are coupled together via (i) nonlinear coupling with coefficients 
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g n _k and (ii) linear coupling with coefficients K n y, where, by symmetry, g n k — g kn 
and K n _k — Hk,n- The nonlinear coupling results from inter-atomic collisions while 
the linear coupling accounts for spin state interconversion usually induced by a spin- 
flipping resonant electromagnetic wave [296]. In the case of fermionic mixtures one 
needs to replace the self-interacting nonlinear terms by <7n.n|V ; n| 4 ^ 3 '0n [297-299]. In 
the absence of losses (cr„ = 0), the total number of atoms is conserved: 



In fact, in the further absence of linear interconversions {n n ^ — 0) each norm N k is 
conserved separately. 

The simplest case of two species (TV = 2) has been studied extensively. In 
particular, if one considers the trapless system (V n = 0) in the absence of linear 
interconversion, losses, and chemical potential differences, the two components tend 
to segregate if the immiscibility condition 



is satisfied [300] . This condition can be interpreted as if the mutual repulsion between 
species is stronger than the combined self-repulsions. In typical experiments, the 
miscibility parameter (an adimensional quantity) is rather small: A m 9 x 10~ 4 for 
87 Rb [289,301] and A w 0.036 for 23 Na [66]. Depending on the various nonlinear 
coefficients, a vast array of solutions can be supported by a binary condensate. 
These include, ground-state solutions [302-304], small-amplitude excitations [305- 
308], bound states of dark-bright [309] and dark-dark [310], dark-gray, bright- 
gray, bright-antidark and dark-antidark [311] complexes of solitary waves, vector 
solitons with embedded domain- walls (DWs) [312], spatially periodic states [313], 
and modulated amplitude waves [314]. Extensions of some of these patterns in 
two dimensions, namely DWs, have been investigated in Refs. [303,304,315-321]. 
The non-equilibrium dynamics of a binary condensate has been shown to support 
(experimentally and theoretically) long lived ring excitations whereby each component 
inter-penetrates the other one repeatedly [322] (see Fig. 8). The effects of adding the 
linear inter-species coupling between the components has also been studied in some 
detail [313,323-333]. One of the salient features of adding the linear inter-species 
coupling is the suppression or promotion of the transition to miscibility (cf. Rcf. [334] 
and Chap. 15 in Ref. [55]). Spinor condensates with three species have also drawn 
considerable attention since their experimental creation [66,335]. Such systems give 
rise to spin domains [290], polarized states [336], spin textures [337], and multi- 
component (vectorial) solitons of bright [338-341], dark [342], gap [343], and bright- 
dark [344] types. 

5.2. Vortices and vortex lattices 

Arguably, one of the most striking nonlinear matter-wave manifestations in BECs is 
the possibility of supporting vortices [345,346] . Vortices are characterized by their non- 
zero topological charge S whereby the phase of the wavefunction has a phase jump 
of 27t5' along a closed contour surrounding the core of the vortex (cf. phase profile 
for a singly charged vortex, S = +1, in the right panel of Fig. 9). Historically, the 
first observation of vortices in BECs was achieved [28] by phase imprinting between 
two hyperfine spin states of Rb [347]. Nowadays, the standard technique to nucleate 




(151) 



A = (.912321 - 911922)/ 9n > 



(152) 
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Figure 8. (Color online). Nonlinear dynamics of a binary 50:50 mixture 
of two spin states (|1, — 1) and 2, 1>) of N = 375,000 Rb atoms. Each 
component is depicted by a contour slice at half of its corresponding maximal 
density. The bottom (side) projection corresponds to the z- (x-) integrated 
density for the |1, — 1) component as it is observed in the laboratory experiments 
[322], Please visit: http://www.rohan.sdsu.edu/~rcarrete/ [Publications] to view 
movies showing the inter-penetrating time evolution between the two components 
over a span of 220 ms. Reprinted from Ref. [322] with permission. 

vortices in BECs is based on stirring [35] the condensate cloud above a certain critical 
angular speed [348-350]. This technique has proven to be extremely efficient, not only 
in creating single vortices, but also, from a few vortices [350], to vortex lattices [351]. 
It is also possible to nucleate vortices by dragging a moving impurity through the 
condensate for speeds above a critical velocity (depending on the local density and 
also the shape of the impurity) [352-359]. Yet another possibility to nucleate vortices 
can be achieved by separating the condensate in different fragments and allowing them 
to collide [360-362]. 




Figure 9. (Color online) Two-dimensional, singly charged (S = 1), vortex inside 
a parabolic magnetic trap V(z) = ^Q 2 (x 2 + j/ 2 ) with f2 = 0.2. Depicted are the 
density (left) and the phase profile (right) for a chemical potential fi = 2. 

The profile of a vortex in a two dimensional setting (see left panel of Fig. 9) can 
be obtained by solving for the density U (r) when considering a wave function of the 
form ?/>(?", 0) = U(r) exp(iS8 — ifit) that satisfies the 2D GP equation with repulsive 
nonlinearity (g = +1), where (r, 9) are the polar coordinates and fi is the chemical 
potential. The equation for U takes the form 

-pi + -^r - V + (" - v ^ - u2 ) u = °> ( 153 ) 

ar z r ar r z 

with boundary conditions U(0) = and U(+oo) = y/Jl for a confining potential 
V(r — +oo) = +oo. Unfortunately, the ensuing equation for the vortex profile cannot 
be solved exactly (even in the simplest homogeneous case with V(r) = 0) and one 
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Figure 10. (Color online). Top: Experimental look at component separation in 
a rotating BEC of Rb atoms [369] after 30 ms evolution, release and 20 ms free 
expansion from a relatively tight trapping potential. The percentages quoted are 
the fraction of atoms transferred from the |1, —1) state (top row) to the [2, 1) state 
(bottom row). Bottom: Numerical vortex lattice (VL) transfer by linear coupling 
from the first (top row) to the second component (bottom row). The initial VL 
(left column) is successfully transferred between components (see middle column 
where the scale in the top panel clearly indicates that the first component is 
almost depleted of atoms after transfer). More importantly, note that the phase 
distribution is also transferred between components (right column). 

has to resort to numerical or approximate methods [363] . The asymptotic behavior of 
the vortex profile U(r) can be found from Eq. (153), i.e., U(r) ~ r' s l as r — ► 0, and 
U(r) ~ JJi — S 2 / (2r 2 ) as r — * +oo. Note that the width of singly charged vortices in 
BECs is of 0(£) [where £ is the healing length given in Eq. (19)], while higher-charge 
vortices (\S\ > 1) have cores wider than the healing length and are unstable in the 
homogeneous background case (V(r) = 0) but might be rendered stable by external 
impurities [364] or by external potentials [365-368]. When unstable, higher order 
charge vortices typically split in multiple single charge vortices. 

Single charge vortices are extremely robust due to their inherent topological 
charge since continuous transformations/deformation of the vortex profile cannot 
eliminate the 2irS phase jump — unless that the density is close to zero (this is the 
reason why, in the stirring experiments, vortices are nucleated at the periphery of the 
condensate cloud where the density tends to zero for confining potentials [370-373]). 
Vortices are prone to motion induced by gradients in both density and phase of the 
background [374]. These gradients can be induced by an external potential or the 
presence of another vortex. The effect of vortex precession induced by the external 
trap has been extensively studied [375-382]. The motion induced on a vortex by 
another vortex is equivalent to the one observed in fluid vortices whereby vortices 
with same charge travel parallel to each other at constant speed, while vortices of 
opposite charges rotate about each other at constant angular speed. 

Another topic that has attracted an enormous deal of attention in recent years 
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is the ubiquitous existence of robust vortex lattices in rapidly rotating condensates 
consisting of ordered lattices of vortices arranged in triangular configurations (the so- 
called Abrikosov lattices [383] ) . The first experimental observation of vortex lattices 
consisted of just a few (< 15) vortices [351] but soon experiments were able to maintain 
vortex lattices with some 100 vortices [384]. Alternative methods to describe vortex 
lattice configurations have been given in terms of Kelvin's variational principle [385] 
and through a linear algebra formulation [386]. Another approach is to treat each 
vortex as a quasi-particle and apply ideas borrowed from molecular dynamics to 
find the most energetically favorable configurations [387]. Some other interesting 
phenomenology of vortex lattices includes the excitation of Tkachenko modes [388] 
via the annihilation of a central chunk of vortex lattice matter through a localized 
laser heating [389]. 

Yet another promising avenue of research that is currently being explored is the 
topic of vortex lattices in multicomponcnt condensates. For example, starting with 
a two-component BEC mixture with only one atomic species containing a vortex 
lattice, and subsequently "activating" the linear coupling, it is possible to entirely 
transfer the vortex lattice to the second component (cf. results in Fig. 10). This 
"Rabi oscillation" between atomic species [390,313] is an extremely useful tool for 
controllably transferring desirable fractions of atoms from one state to another and can 
be extended to multicomponcnt, condensates [391,392]. Furthermore, it is important 
to note that the interaction of vortex lattices in a multicomponcnt BEC can result in 
structural changes in the configurations of the vortex lattices, i.e., resulting in lattices 
with different symmetry [393,394]. 

5.3. Shock waves 

One of the classical types of nonlinear waves appearing in the context of BECs is shock 
waves. Shock waves were first observed in the experiments reported in Ref. [31], where 
a slow- light technique was used to produce density depressions in a sodium BEC. More 
recently, they were observed in rapidly rotating 87 Rb BECs triggered by repulsive 
laser pulses [395] , while their formation was discussed in an experiment involving the 
growth dynamics of a ID sodium quasi-condensate in a dimple microtrap created on 
top of the harmonic confinement of an atom chip [396]. Finally, shock waves were 
studied experimentally and theoretically in Ref. [397]; in the experiments reported in 
this work, repulsive laser beams were used (as in Ref. [395], but with a nonrotating 
condensate) to push atoms from the BEC center, thus forming "blast-wave" patterns. 

On the theoretical side, shock waves in repulsive BECs were mainly studied in the 
framework of mean-field theory and the GP equation for weakly-interacting Bose gases 
[397-405] , but also for strongly interacting ones [406] and in the BEC-Tonks crossover 
[407]; additionally, the effect of temperature (see Sec. 5.7) on shock wave formation and 
dynamics and the effect of depleted atoms were respectively discussed in Refs. [396] 
and [407] . Many of the above mentioned theoretical studies rely on the hydrodynamic 
equations that can be obtained from the GP equation via the Madelung transformation 
i\) — y/p exp(i</>) (with p and cf) being the condensate's density and phase, respectively, 
see also Sec. 4.3.3). These hydrodynamic equations are then treated in the long- 
wavelength limit (or, equivalently, in the weakly dispersive regime) where the so-called 
quantum pressure term, ~ (V 2 ^)/^, is negligible. Qualitatively speaking, one may 
ignore this term if the so-called "quantum Reynolds number" [399,400] (see also the 
discussion in Refs. [398,402]) is R = anoL/Q 3> 1, where a is the scattering length, 
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n the peak density of the BEC and L a minimal scale among all the characteristic 
spatial scales of the condensate wavcfunction. In fact, a more rigorous treatment 
relies on the assumption that the quantum pressure term is small, as in the case of the 
theoretical analysis (and the pertinent experimental results) of Ref . [397] . In any case, 
since the quantum pressure is reminiscent of viscosity in classical fluid mechanics, this 
weakly dispersive regime suggests the possibility of "dissipationless shock waves" in 
BECs (according to the nomenclature of Ref. [403]). 

The above mentioned hydrodynamic equations were treated in the limiting case 
of zero quantum pressure in Refs. [398-401]. In this case, the hydrodynamic equations 
resulting from the GP equation are reduced to a hyperbolic system of conservation 
laws of classical gas dynamics, namely an Euler and a continuity equation (see, 
e.g., Ref. [408]). In this gas dynamics picture, the above hyperbolic system is 
characterized by two real eigenspeeds (i.e., characteristic speeds of propagation of weak 
discontinuities) = vie, where v is the velocity of the gas and c is the speed of sound. 
Since the latter depends on the condensate density p [see Eq. (18)], it is clear that 
higher values of p propagate faster than lower ones and, as a result, any compressive 
part of the wave ultimately breaks to give a triple-valued solution for p; this is a 
signature of the formation of a shock wave. In this system, a Gaussian input produces 
two symmetric shocks at finite time, as was shown e.g., in Refs. [401,407] (see also 
Ref. [409] for a rigorous discussion) . Importantly, the trailing edge of the shock wave as 
observed in the simulations and the experiments (see, e.g., Refs. [31,396,397,401,407]) 
can be considered as a modulated train of dark solitons [397,403]. 

In the work [402], a more careful investigation of the regimes of validity of the 
condition R » 1 (which may fail, e.g., in the case of expanding BECs) led to an 
experimentally relevant protocol to produce shock waves in BECs. This protocol is 
based on the use of Feshbach resonance to control the scattering length a, namely 
to make a an increasing function of time (by a proper ramp- up procedure), so as 
to increase the time domain in which the quantum pressure is negligible. This way, 
this "Feshbach resonance management" technique ff was shown to produce multi- 
dimensional shock waves. On the other hand, in Ref. [403] the Whitham averaging 
method [413] was used in the weakly dispersive regime to show that dissipationless 
shock waves are emanating from density humps in repulsive BECs. Moreover, the 
formation of shock waves in BECs confined in optical lattices was discussed in 
Ref. [404], while in Ref. [397] an in depth analysis of the shock waves appearing 
in BECs and in gas dynamics (also in connection to relevant experiments of the JILA 
group) was presented. Finally, it should also be mentioned that the above mentioned 
works chiefly refer to repulsive BECs; an analysis of the shock wave formation in 
attractive BECs can be found in Ref. [91]. 

5.4- Multidimensional solitons and collapse 

As was highlighted in Sec. 3, a strong transverse confinement may effectively render 
the BEC quasi-lD, in which case the ID soliton solutions are physically relevant 
and arc stable. On the contrary, in the absence of a tight transversal trapping, 
higher dimensional extensions of ID solitons are generally unstable [12]. However, by 
restricting the transverse direction(s) of the condensate, it is possible to obtain higher 

ff This technique was suggested in earlier works in various important applications, as e.g., a means 
to prevent collapse of higher-dimensional attractive BECs [193,203,410], to produce periodic waves 
[107], robust matter-wave breathers [241,242,411,412], and so on. 
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dimensional soliton solutions that are stabilized for times longer than the lifespan of 
the experiments [414-416]. 

Let us now showcase some of the possible higher dimensional soliton solutions 
displayed by the GP model. We do not cover here "true" 3D solutions (i.e., solutions 
that do not have a ID equivalent) such as 2D vortices (see previous section), 3D vortex 
lines [11,417-420], vortex rings [421,422,33] or more complicated topologically charged 
structures such a skyrmions [423-428]. 

5.4-1- Dark solitons. The trademark of a dark soliton is its phase jump along its 
center separating two repelling phases. In a 2D geometry a dark soliton corresponds 
to a nodal line separating the two phases while in 3D it corresponds to a nodal plane. 
Both the 2D and 3D dark solitons, respectively called band (or stripe) and planar 
dark solitons, are prone to the snaking instability along their nodal extent [429,430]. 
These instabilities result in the nucleation of vortex pairs in 2D [431,432] and pairs of 
vortex lines and/or vortex rings in 3D. When the dark soliton is set into motion inside 
a confining trap, it suffers bending resulting from the different speeds of sounds at the 
edge of the cloud (low density and thus slower speeds compared to the speed at the 
center of the cloud) accelerating the formation of vortices at the trailing edges [30] . It 
is also possible to create dark soliton structures whose nodal sets, instead of extending 
linearly, can be wrapped around. It is therefore possible to create in 2D ring dark 
solitons and in 3D spherical shell dark solitons, as the ones discussed in Sec. 4.3.2. 
Such structures can be described by nonlinear Bessel functions (cf. Ref. [433] and 
Chap. 7 in Ref. [55]), are also prone to the snaking instability. It is worth mentioning 
that the abovementioned instabilities can be weak (slow) enough so that these solitons 
can be observed in the experiments [33,31,434]. 

5-4-2. Bright solitons and collapse. Bright solitons in higher dimensions are prone 
to a different type of instability due to the intrinsic collapse of solutions of the NLS 
equation [12]. The first experiments with attractive condensates suffered from this 
collapse instability [81,82,435,436] while the more recent experiments were able to 
focus on stable regimes and demonstrate bright soliton formation [40-42]. In fact, 
the key feature of these experiments was the quasi-lD nature of the attractive BEC 
realized in anisotropic traps as discussed in Sec. 3. Thus, the observed bright solitons 
were found to be robust, which would not be the case in a higher-dimensional system, 
as they should either collapse or expand indefinitely depending on the number of atoms 
and the density profile. The solution that constitutes the unstable separatrix between 
expansion and collapse is the well-known Townes soliton [437,438]. 

In this connection, it is important to mention that even though the experimental 
condensates are never truly ID, fortunately, the tight trapping in the transverse 
direction(s) is able to slow collapse to times much longer than the duration of the 
experiments. Nonetheless, interactions of bright solitons in higher dimensions, in 
contrast with their ID counterparts, may be inelastic [117,89,439], and, furthermore, 
when two (or more) solitons overlap their combined number of atoms can exceed the 
critical threshold and initiate collapse. The overlap of higher dimensional solitons, 
even when the critical number of atoms is exceeded, might not result in collapse since 
one has to take into account the time of the interaction (depending on the velocities of 
the solitons and their relatives phases, cf. Chap. 7 in Ref. [55] and references therein). 

Finally, we note that stabilization of higher-dimensional bright solitons by 
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means of lower-dimensional optical lattices has been proposed in Refs. [192,440,412]. 
Moreover, stable bright ring solitons carrying topological charge have been 
theoretically predicted to exist [441,83,442,443,366] but no experiment has yet 
corroborated these results. 

5.5. Manipulation of matter-waves 

As discussed in Sec. 3, one of the most appealing features in BECs is the level of control 
over the different contributions on the GP equation. This includes the possibility 
to craft almost any desired external potential (by the appropriate superposition of 
multiple laser beams), and to change the strength and sign of the nonlinearity via 
the Feshbach resonance mechanism. This is to be contrasted with other contexts 
where the NLS equation is also a relevant model, as, e.g., in nonlinear optics [17], 
where it is extremely difficult and often impossible to demonstrate such control. 
In the BEC context, particularly appealing is the fact that the external potential 
and/or nonlinearity can be made to follow in time any desired evolution. In this 
section we focus on the use of appropriately crafted time-dependent external potentials 
to manipulate mater-waves in BECs. In this section we only consider matter- 
wave manipulation by two main types of external potentials: localized and periodic 
potentials. These potentials are experimentally generated by laser beams as explained 
in Sec. 2.4. Here, it would be useful to recall that the sign of the external optical 
potential is positive (negative) for blue- (red-) detuned laser beams. 

5.5.1. Localized potentials. The interaction of solitons with localized impurities has 
attracted much attention in the theory of nonlinear waves [253,254,444] and solid 
state physics [445,446]. In ID BECs confined by a harmonic trap, bright and dark 
solitons perform harmonic oscillations as discussed in Sec. 4 as a consequence of 
the Kohn's theorem (the "nonlinear analogue" of the Ehrcnfest theorem) [198,447], 
which states that the motion of the center of mass of a cloud of particles trapped 
in a parabolic potential decouples from the internal excitations. The existence, 
stability and dynamics of bright solitons in the presence of the external potential 
can be analyzed using perturbation techniques expounded in Sec. 4. In fact, the 
combined effects of the harmonic trap Vht(^) = Sl 2 x 2 /2 and an infinitely localized 
delta impurity, located at Xi, namely Vi mp (x) = 8{x — Xi), yield the following 
effective force on a bright soliton [448], 

F cS = F HT + Fimp - -2ft 2 ?? ( - 2r] 3 V^l tanh(r/(a; 4 - 0) scch 2 ^ - C)), (154) 

where £ and rj are, respectively, the center and height of the bright soliton, while 
the first (second) term in the right-hand side corresponds to the harmonic trapping 
(localized impurity) potential. This effective force induces a Newton-type dynamics 
for the center £ of the bright soliton, as discussed previously. Note that the force 
induced by the harmonic trap always points towards its center, while the direction 
of the force induced by the impurity depends on the sign of V^ p . In the case of an 
attractive impurity it is possible to not only pin the bright soliton away from the center 
of the harmonic trap but also to adiabatically drag it and reposition it at, almost, any 
desired location by slowly moving the impurity [448]. The success in dragging the 
soliton not only depends on the profiles of the soliton and the impurity, but more 
crucially, on the degree of adiabaticity when displacing the impurity. 
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A similar study can also be performed in the case of dark solitons. The 
interactions of dark solitons with localized impurities was analyzed in Ref. [449], by 
using the so-called direct perturbation theory for dark solitons [450], and later in the 
BEC context in Ref. [258] , by using the adiabatic perturbation theory for dark solitons 
[114]. Following the analysis of Ref. [258], it is possible to show that the center £ of 
a dark matter- wave soliton obeys a Newton-type equation of motion with an effective 
potential given by: 

Kff(C) = V3?(C) + ^imp(C) = ^Vht(C) + J^sech 2 (C). (155) 

Thus, similarly to the bright soliton case, one can use the above pinning of the dark 
soliton to drag it by adiabatically moving the impurity [258] . 

Finally, it is also possible to pin and drag vortices with a localized impurity 
(see Chap. 17 in Ref. [55] for more details). As discussed above, the presence of the 
harmonic trap induces vortex precession [375-382]. Therefore, in order to pin/drag 
the vortex at an off-center position, the pinning force exerted by the impurity has to 
be stronger than the vortex precession force induced by the harmonic trap and the 
impurity has to be deep enough to avoid emission of sound waves [451]. An interesting 
twist to this manipulation problem is the possibility to snare a moving vortex by an 
appropriately located/designed impurity. 

5.5.2. Periodic potentials. A similar approach as the one described in the previous 
section can be devised by using periodic potentials. This method relies, as in the 
case of localized impurities, on the pinning properties of properly crafted periodic 
potentials. Specifically, a periodic potential with a wavelength longer than the width 
of the soliton width induces a periodic series of effective potential minima that help 
to pin the soliton. For instance, in a ID BEC, a bright soliton subject to a periodic 
potential generated by an optical lattice of the form 

Vol{x) =V^sm 2 {kx), (156) 

behaves (in appropriate parameter ranges) as a quasi-particle inside the effective 
potential [257]: 

Kff(C) = V&C 2 - 7rV^L ) fccsch(fc7r/77)cos(2fcC). (157) 

This effective potential possesses minima that, as indicated above, can be rigorously 
shown to correspond to stable positions for the solitons [217-222,452,453] . This in turn 
can be used, in the same manner as above for the localized impurities, to successfully 
pin, drag and capture bright solitons [257]. 

The case of manipulation of dark solitons by periodic potentials is more subtle 
because dark solitons subject to tight confinements are prone to weak radiation loss, 
as shown numerically in Refs. [454-457] and analytically in Ref. [263]. Since a dark 
soliton is an effectively negative mass (density void) structure, radiation loss implies 
that the soliton travels faster and eventually escapes the local minimum in the effective 
potential landscape. The motion of a dark soliton subject to an external potential has 
been treated in detail before [258-263,458]. In the presence of both the magnetic trap 
and the optical lattice of Eq. (156), the motion of the dark soliton depends on the 
period of the optical lattice when compared to the soliton's width [459] . The case of 
an optical lattice with long-period can be treated as a perturbation (see Sec. 4.3.2), 
and the dynamics of the dark soliton can adjust itself to the smooth potential. The 
short period optical lattice case can be treated by multiple-scale expansion [459] and 
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it is equivalent to the motion of a dark soliton inside a renormalized magnetic trap 
(with no optical lattice) [459]. For intermediate periods, the optical lattice has the 
ability to drag/manipulate the dark soliton as shown in Rcfs. [459,460]. However, the 
effects of radiation loss described above eventually drive the dark soliton into large 
amplitude oscillations inside the local effective potential minima leading, eventually, 
to its expulsion. 

Finally, we briefly discuss vortices under the presence of periodic potentials. Their 
stability and dynamics has been studied in Refs. [461,462], while the existence of gap 
vortices in the gaps of the band-gap spectrum due to Bragg scattering were considered 
in Ref. [177]. Also, in the presence of deep periodic lattices, where the discrete version 
of the GP equation (namely the DNLS equation, see, e.g., Sec. 3.6.2) becomes relevant, 
purely 3D discrete vortices can be constructed [463-465]. An interesting twist of the 
pinning of vortices (which has been observed experimentally [466]) is the case when 
a vortex lattice is induced to transition from a triangular Abrikosov vortex lattice 
to a square lattice by an optical lattice [467,468]. Another type of vortex lattice 
manipulation is the use of large-amplitude oscillations to induce structural phase 
transitions (e.g., from triangular to orthorhombic) [39]. 

5.6. Matter-waves in disordered potentials 

Recently, there has been much attention focused on the dynamics of matter-waves in 
disordered potentials. Generally speaking, disorder in quantum systems has been a 
subject of intense theoretical and experimental studies. In the context of ultracold 
atomic gases disorder may result from the roughness of a magnetic trap [469] or a 
magnetic microtrap [470]. However, it is important to note that controlled disorder (or 
quasi-disorder) may also be created by means of different techniques. These include 
the use of two-color supcrlattice potentials [471-473], the employment of so-called 
quasi-crystal (i.e., quasi-periodic) optical lattices in 2D or 3D [474-476], the use of 
impurity atoms trapped at random positions in the nodes of a periodic optical lattice 
[477], random phase masks [478], or optical speckle patterns [479-481]. The latter is 
a random intensity pattern which is produced by the scattering of a coherent laser 
beam from a rough surface (see, e.g., Ref. [482] for a detailed discussion). 

The theoretical and experimental investigations of ultracold atomic gases confined 
in disordered potentials pave the way for the study of fundamental effects in quantum 
systems. Among them, the most famous phenomenon is Anderson localization, 
i.e., localization and absence of diffusion of non- interacting quantum particles, 
which was originally predicted to occur in the context of electronic transport [483]. 
Other important effects, include the realization of Bose [484,485] or Fermi [486,487] 
glasses, quantum spin glasses [488], the Anderson-Bose glass and crossover between 
Anderson glass to Bose-glass localization, and so on (see also the recent review [489]). 
Importantly, as we will discuss below, there exist very recent relevant experimental 
results (and theoretical predictions) towards these interesting directions. 

The first experimental results concerning BECs confined in disordered potentials 
were reported almost simultaneously on 2005 [478-481]. In the first if these 
experiments [479] , static and dynamic properties of a rubidium BEC were studied and 
it was found that both dipole and quadrupole oscillations are damped (the damping 
was found to be stronger as the speckle height was increased). The suppression of 
transport of the condensate in the presence of a random potential was also reported 
in Refs. [480,481]. In these works, a strong reduction of the mobility of atoms was 
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demonstrated, as the ID expansion of the elongated condensate along a magnetic 
[480] or an optical [481] guide was found to be strongly inhibited in the presence of the 
speckle potential. On the other hand, in Ref. [478] a speckle pattern was superimposed 
to a regular ID optical lattice and, thus, a genuine random potential was created. In 
this setting the possibility of observation of Anderson localization was analyzed in 
detail, and the crossover from Anderson localization in the absence of interactions to 
the "screening regime" (where nonlinear interactions suppress Anderson localization 
in the random potential) was investigated. 

Although the above mentioned results underscore a disorder-induced trapping of 
the condensate, this effect does not correspond to a genuine Anderson localization: 
for the latter, the correlation length of the disorder has to be smaller than the size 
of the system (see discussion in Ref. [489]). Nevertheless, a detailed study [490] 
has shown that expansion of a quasi-2D cloud may lead to weak and even strong 
localization using currently available speckle patterns. Other works have revealed that 
Anderson localization may occur during transport processes in repulsive BECs [491- 
493]. In this case, however, and for condensates at equilibrium, the interaction-induced 
derealization dominates disorder-induced localization, except for the case of weak 
interactions [494] (see also earlier work in Ref. [495]). Another possibility is Anderson 
localization of elementary excitations in interacting BECs, as analyzed in the recent 
works [496,497]. Finally, it is worth also mentioning in passing parallel developments 
in this area, within the mathematically similar setting of photonic lattices, where 
Anderson localization and transition from ballistic to diffusive transport were recently 
observed in the presence of random fluctuations [498] . 

In any case, the above discussion shows that there exists an intense theoretical 
and experimental effort concerning this hot topic of BECs in disordered potentials. 
Although it seems that relevant experiments have just started, the perspectives are 
very promising: they include not only the possibility of the detection of Anderson 
localization, but also other relevant effects, such as the observation of a Bose-glass 
phase [499] , the possibility of the appearance of a novel Lifshits glass phase [500] , and 
so on. 

5.7. Beyond mean-field description 

The GP equation has been extremely successful in describing a wide range of mean- 
field phenomena in Bose-Einstein condensation. By construction, as explained in detail 
in Sec. 2, the GP equation is the mean-field description of the multi-body quantum 
Hamiltonian describing the interaction of a dilute gas of bosonic atoms and it relies on 
two main assumptions: (a) collisions between atoms are approximated by hard sphere 
collisions with a Dirac delta interatomic potential, and (b) the gas is at absolute 
zero temperature where thermal effects are not present. Nonetheless, in many BEC 
settings, finite temperature effects and quantum fluctuations may play an important 
role. The main effect of finite temperature is due to the fact that a part of the atomic 
cloud is not condensed (the so-called thermal cloud) and couples to the condensed 
cloud. A microscopic derivation of the mean-field operator for the gas of bosonic 
particles reveals that the (standard) condensate mean-field is coupled to higher order 
mean-fields (cf. the insightful review in Chap. 18 of Ref. [55] and references therein 
for more details) . Neglecting the exchange of particles between condensed and non- 
condensed atoms and taking into account the three lowest mean- field orders, the so- 
called Hartree-Fock-Bogoliubov (HFB) theory, leads to the generalized GP equation 
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for the condensate wavefunction tl>(r,t) [501-504] 

ih^{r,t) = [H GP + g[2n'(r,t)}}iP(r,t) +m (r,t)i>*{r,t), (158) 

where H GP = -(1/2)V 2 + V cx t(Y,t) + g\ip\ 2 accounts for the "classical" GP terms 
in non-dimensional units, n'(r,t) denotes the non-condensate density and m (r, t) 
the anomalous mean-field average [505-508]. Furthermore, taking into account that 
atomic collisions happen within the gas (and not in vacuum), one has to modify the 
inter-atomic interactions by a contact potential with a position-dependent amplitude: 
g -» g[l + m (r)/^ 2 (r)] [509-511]. 

A similar approach to the above is to consider the coupling with the thermal 
cloud by ignoring the anomalous average and including a local energy and momentum 
conservation [512,513] that yields 

ih^(r, t) = [H GP + g [2n'(r, t)] - iR(r, t)} ^(r, t), (159) 

where the non-condensate density n'(r, t) is described terms of a Wigner phase-space 
representation and a generalized Boltzmann quantum-hydrodynamical equation (see 
Chap. 18 of Ref. [55]). This term includes collisions between non-condensed atoms and 
transfer to and from the condensed cloud. In Eq. (159), R{r,t) accounts for triplet 
correlations. 

Other approaches to incorporate the thermal cloud include Stoof's non- 
equilibrium theory based on quantum kinetic theory [514-517] and Gardiner-Zoller's 
quantum kinetic master equation using techniques borrowed from quantum optics 
[518-524]. Also, efforts to include stochastic effects have been considered in Ref. [525] 
by considering the thermal cloud to be close enough to equilibrium that it can be 
described by a Bose cloud with chemical potential fir at temperature T. In fact, 
this approach leads, after neglecting noise terms, to the phenomenological damping 
included in the GP equation through the damped GP equation 

i»J^(r, *) = (!- h) (Hop - Mr) ^(r, t), (160) 

with damping rate 7 > 0. This approach was originally proposed by Pitaevskii 
[526] and applied with a position-dependent loss rate in Ref. [527]. A similar 
phenomenological damping term (but on the left hand side) has been used in 
Refs. [371,372] to remove the excess energy to dynamically simulate the crystallization 
of vortex lattices in rapidly rotating BECs. 

It is important to note that recently there has been a lot of activity on the study 
of formation, stability and dynamics of matter-wave solitons beyond the mean-field 
theory. In particular, in the case of attractive BECs, the velocity and temperature- 
dependent frictional force and diffusion coefficient of a matter-wave bright soliton 
immersed in a thermal cloud was calculated in Ref. [116]. Moreover, the full set of the 
time-dependent HFB equations was used in Ref. [528] to show that the matter flux 
from the condensate to the thermal cloud may cause the bright matter-wave solitons 
to split into two solitonic fragments (each of them is a mixture of the condensed 
and non-condensed particles); these may be viewed as partially incoherent solitons, 
similar to the ones known in the context of nonlinear optics [529,530]. Partially 
incoherent lattice solitons at a finite temperature T were also predicted to exist in 
Ref. [531]; there, a numerical integration of the GP equation, starting with a random 
Bose distribution at finite T, revealed the generation of lattice solitons upon gradual 
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switch-on of an optical lattice potential. In a more recent study [532], the time- 
dependent HFB theory was used to find inter-gap and intra-gap partially incoherent 
solitons (composed, as in Ref. [528] by a condensed and a non-condensed part). On the 
other hand, in Refs. [533,534], the so-called Zarcmba-Nikuni-Griffin formalism [512] 
was used to analyze the dissipative dynamics of dark solitons in the presence of the 
thermal cloud, as observed in the Hannover experiment [32]. Finally, other quantum 
effects associated to matter-wave solitons have been considered as well; these include 
quantum depiction of dark solitons [535-538] (see also Ref. [115]), formation of a 
broken-symmetry bright matter-wave soliton by superpositions of quasi-degenerate 
many-body states [539], quantum-noise squeezing and quantum correlations of gap 
solitons [540], and so on. 

A considerable volume of work in many of the directions mentioned in this review 
has continued to emerge both in preprint and in published form. As two examples 
thereof, we can mention the work of Refs. [541-543] on various higher dimensional 
aspects of BECs in random potentials and transitions associated with Anderson 
localization (we thank Dr. B. Shapiro for bringing these works to our attention) and 
the intense experimental activity on the oscillations and interactions of dark and dark- 
bright solitons in Refs. [544,545]. 
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